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the  present  findings,  we  present  an  algorithm  which  permits  recovery  of  the  steady  state  from 
transient  data  at  a  negligible  fraction  of  the  otherwise  necessary  computing  effort.  Application 
of  the  algorithm  to  extract  the  three-dimensional  flowfield  on  the  elliptic  cone  is  underway. 
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1.  Introduction 

Revived  interest  in  affordable  hypersonic  flow  technology  (Walker  1999)  has  resulted  in  re¬ 
newed  efforts  to  understand  and  control  the  related  physical  mechanisms  of  laminar-turbulent 
flow  transition  (Kimmel  and  Walker  1999).  One  of  the  configurations  in  discussion  is  that  of 
the  elliptic  cone  which  is  known  to  have  distinct  aerodynamic  advantages  over  bodies  of  revolu¬ 
tion  in  compressible  flow  at  least  as  early  as  the  pioneering  work  of  Jorgensen  (1958).  Renewed 
theoretical  and  experimental  efforts  on  elliptic  cones  have  been  undertaken  by  Kimmel  et  al. 
(1997,1999)  and  Huntley  et  al.  (1999)  with  special  emphasis  placed  on  laminar-turbulent  flow 
transition,  thus  expanding  upon  knowledge  on  compressible  boundary-layer  flow  transition  on  a 
circular-base  cone  (Stetson  and  Kimmel  1992).  Kimmel  et  al.  (1999)  have  discovered  character¬ 
istics  of  flow  transition  on  the  elliptic  cone  which  appear  to  be  intrinsically  three-dimensional 
and  distinct  from  known  instability  mechanisms  on  bodies  of  revolution.  Despite  substantial 
progress  on  the  issue,  three-dimensional  boundary  layer  transition  in  high-speed  flow  presently 
appears  to  be  far  from  understood  in  a  satisfactory  manner. 

The  issue  of  laminar-turbulent  transition  has  a  long  and  successful  track  record,  following  the 
pioneering  works  of  the  first  half  of  last  century  on  linear  instability  in  boundary  layers  on  flat- 
plates  and  bodies  of  revolution  (Lin  1955).  The  main  limitation  of  most  of  the  existing  analyses  is 
their  confinement  within  the  framework  of  flowfields  in  which  two  of  the  three  spatial  directions 
are  taken  to  be  homogeneous  and  are  resolved  by  a  periodicity  Ansatz.  Under  these  so-called 
local  assumptions  the  classic  Rayleigh  and  Orr-Sommerfeld  equations  result,  respectively,  in 
inviscid  and  viscous  incompressible  flow  (Drazin  and  Reid  1981).  The  compressible  counterparts 
of  either  equation  have  been  extensively  investigated  and  discussed  by  Mack  (1969,  1984).  With 
the  maturing  of  numerical  algorithms  and  the  current  progress  in  hardware  technology  the 
assumption  of  homogeneity  in  two  spatial  directions  can  be  relaxed  within  the  framework  of 
a  direct  numerical  simulation  (DNS)  in  which  all  three  spatial  directions  may,  in  principle,  be 
taken  as  inhomogeneous  and  be  resolved;  marching  the  equations  of  motion  forward  in  time 
delivers  the  desired  instability  information.  While  DNS  is  an  undisputed  approach  in  terms  of 
its  capacity  to  deliver  solutions  of  the  equations  of  motion,  there  are  two  points  that  deserve 
discussion  with  respect  to  using  DNS  for  instability  analyses. 

The  first  point  is  related  with  computing  effort  and  the  associated  cost  of  a  DNS.  The  latter, 
though  decreasing  in  recent  years  in  proportion  to  the  continuous  increase  in  hardware  capabili¬ 
ties,  is  still  formidable;  the  storage  requirements  for  resolution  of  all  scales  of  three-dimensional 
turbulent  flow  scale  with  the  9/4  —  th  power  of  an  appropriately  defined  flow  Reynolds  num¬ 
ber.  As  a  matter  of  fact,  the  resolution  requirements  for  a  simulation  through  laminar-turbulent 
transition  into  turbulent  flow  are  still  higher  (e.g.  Gilbert  and  Kleiser  1988)  owing  to  the  steep 
gradients  developing  during  the  multi-spike  stages  of  transition.  Furthermore,  compressibility 
generally  lowers  the  growth  rates  of  the  most  amplified  linear  disturbances  in  comparison  with 
those  in  incompressible  flow.  Consequently,  instability  analyses  of  three-dimensional  compress¬ 
ible  flows  using  DNS  are  still  rare  and  when  performed  are  restricted  to  a  small  subspace  of  pa¬ 
rameters.  Notable  examples  of  compressible  flow  instability  analyses  based  on  DNS  have  started 
to  appear  in  the  literature  (e.g.  Pruett  and  Zang  1992,  Zhong  1998).  Their  obvious  strength 
is  that  both  linear  and  nonlinear  stages  of  flow  instability  can  be  tackled  in  a  self-consistent 
manner. 

The  second  point  is  rather  more  subtle  and  is  related  with  the  potential  of  confusion  of  the 
different  types  of  local  or  global  linear  instabilities,  which  potentially  coexist  in  the  problem 
and  are  all  resolved  by  the  DNS.  The  self-contained  work  to  be  found  as  section  6  of  this 
report  elucidates  this  point  in  detail.  Here  it  suffices  to  stress  the  difference  between  an  ordinary 
differential  equation  based  linear  theory,  such  as  an  Orr-Sommerfeld  type  of  analysis,  and  one  in 
which  the  partial  differential  equation  eigenvalue  problem  is  solved.  Another  point  that  deserves 
mention  here  is  that  for  flows  with  a  mild  variation  in  one  spatial  direction  Herbert  (1997  and 
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references  therein)  has  proposed  the  concept  of  parabolised  stability  equations  (PSE)  as  an 
alternative  tool  to  both  an  Orr-Sommerfeld  type  of  analysis  and  DNS.  Using  the  PSE  linear 
and  nonlinear  instability  analysis  results  inaccessible  to  the  classic  linear  theory  (and  obtainable 
by  DNS)  can  be  successfully  recovered  at  a  negligible  fraction  of  the  effort  required  by  DNS. 
However,  when  the  basic  flow  to  be  analysed  is  essentially  two-  or  three-dimensional  PSE  is 
inapplicable  to  its  instability  analysis.  There  is  clearly  room  for  instability  analysis  tools  beyond 
a  local  linear  theory  based  on  an  Orr-Sommerfeld  type  of  equation  and  alternative  to  both  DNS 
and  PSE. 

Such  an  analysis  tool  may  be  developed  within  the  framework  of  a  global  linear  instability 
theory  the  principles  of  which  are  discussed  in  section  4  and  detailed  in  section  6  b  (iii).  This 
theory  is  analogous  in  spirit,  if  not  in  detail,  to  the  local  linear  analysis  in  considering  the  flowfield 
as  composed  of  a  steady  (or  time-periodic)  so-called  ’basic-flow’  part  upon  which  small-amplitude 
three-dimensional  disturbances  are  superimposed.  The  key  difference  between  global  and  local 
theory  is  that  in  the  former  both  the  basic  flow  and  the  amplitude  functions  of  the  linear 
perturbations  are  essentially  two-  or  three-dimensional  nonperiodic  functions  of  the  respective 
resolved  spatial  coordinates.  The  resulting  freedom  with  the  boundary  conditions  to  be  imposed 
on  the  elliptic  global  problem  permits  the  recovery  of  much  wider  classes  of  linear  instabilities 
when  compared  with  those  delivered  by  the  local  analysis.  The  development  of  the  global  linear 
disturbances  in  time  is  studied  through  the  solution  of  a  large  non-symmetric  complex  generalised 
eigenvalue  problem  in  which  the  governing  equations  may  be  recast  upon  substitution  of  the 
decomposition  and  linearisation  about  the  basic  state.  Owing  to  the  orders-of-magnitude  higher 
computing  effort  compared  with  the  classic  linear  analysis,  global  linear  analysis  results  are 
slowly  emerging  in  the  literature,  mainly  confined  to  flows  possessing  special  symmetries  which 
reduce  the  computational  effort  related  with  the  solution  of  the  eigenvalue  problem  (e.g.  Tatsumi 
and  Yoshimura  1990,  Hall  and  Horseman  1991,  Lin  and  Malik  1996).  An  incompressible  global 
linear  analysis  algorithm  was  developed,  in  primitive  variables  and  without  the  restriction  to 
flows  with  special  symmetries,  by  Theofilis  (1998c);  the  algorithm  has  been  employed,  amongst 
other  flows,  to  the  discovery  of  global  instability  in  separated  flat-plate  boundary  layer  flow 
(Theofilis  et  al.  2000)  while  its  predictions  in  the  classic  lid-driven  cavity  (Theofilis  2000)  are  in 
excellent  agreement  with  the  experiments  of  Aidun  et  al.  (1991)  and  Benson  and  Aidun  (1992). 
A  recent  review  of  the  global  linear  analysis  results  is  presented  by  Theofilis  (2001). 

The  present  report  is  the  first  step  towards  global  linear  instability  analysis  of  three-dimensional 
flowfields  over  elliptic  cones  at  hypersonic  speeds.  To  this  end,  we  discuss  here  the  basic  flows 
which  form  the  variable  coefficients  of  the  partial  differential  equation  eigenvalue  problem  at  a 
small  representative  class  of  the  parameter  values  at  which  the  analyses  will  ultimately  be  per¬ 
formed.  The  essentials  of  the  numerical  issues  surrounding  the  generation  of  the  basic  flows  are 
discussed  first  and  laminar  steady-state  solutions  of  the  equations  of  motion  are  subsequently 
presented.  We  use  this  report  to  illustrate  characteristic  compressible  flowfields  recovered  at  con¬ 
stant  Reynolds  number  and  increasing  angle  of  attack  and  Mach  number,  keeping  for  the  most 
part  of  this  work  the  same  resolution  of  three-dimensional  space.  Grid  adaptation  was  found  to 
be  necessary  as  the  Mach  number  increases,  on  account  of  increasingly  stronger  discontinuities 
and  steeper  gradients  in  the  flow.  The  necessary  but  computationally  intensive  grid  refinement 
studies,  which  will  assess  the  adequacy  of  the  recovered  solutions  for  an  instability  analysisf,  are 
left  to  be  addressed  when  the  analysis  will  be  performed  at  specific  sets  of  parameters. 

The  analytical  description  of  the  elliptic  cone  surface  in  a  body-fitted  cartesian  coordinate 
system  Ox'y'z'  is  provided  by  one  of  the  roots  x'  of  the  equation 


t  f.e.  in  terms  of  convergence  of  second  wall-normal  and  azimuthal  derivatives 
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+  +  ^  =  (1-1) 
a^  0"^  c^ 

where  a,  6,  c  are  the  semiaxes  of  the  elliptic  cone.  The  geometry  of  interest  is  schematically 
depicted  in  figure  1.  In  this  figure  the  angles  of  attack  a  and  bank  are  also  shown.  Interestingly, 
an  effect  of  a  nonzero  angle  on  the  global  instability  analysis  results  is  only  to  be  expected 
when  a  /  0.  In  that  case,  at  each  value  of  the  angle  of  attack  values  of  angle  of  bank  0  <  /3  <  7r/2 
must  be  examined  in  order  for  the  global  instability  on  the  elliptic  cone  to  be  fully  understood.  If 
a  =  0,  on  the  other  hand,  the  basic  flowfields  are  invariant  against  a  change  in  and  the  number 
of  parameters  in  the  problem  is  reduced  by  one.  Another  simplification  of  the  global  eigenvalue 
problem  in  the  limiting  case  a  =  0,  to  be  discussed  in  detail  in  what  follows,  is  the  ability  to 
perform  an  inviscid  global  analysis  in  order  to  extract  the  essential  flow  physics  while  reducing 
the  computing  effort.  This  is  achieved  by  use  of  the  generalised  (within  the  framework  of  global 
linear  theory)  compressible  Rayleigh  equation,  derived  and  presented  for  the  first  time  as  part  of 
the  present  work.  Further,  following  the  recent  discovery  of  the  connection  of  numerical  residuals 
in  steady  flow  calculations  to  global  linear  flow  eigenmodes  it  has  become  clear  that  the  basic 
flow  issue  cannot  be  treated  independently  from  that  of  the  global  linear  instability  analysis,  if 
the  latter  is  to  be  performed  in  a  self-consistent  manner. 

In  the  following  section  2  we  discuss  essential  elements  of  the  approach  taken  for  the  generation 
of  basic  flow  data.  In  section  3  results  on  model  problems,  obtained  by  employing  the  AUSMDV 
scheme  chosen  for  the  generation  of  the  basic  states,  are  presented  before  our  attention  is  focussed 
on  the  three-dimensional  laminar  flowfields  on  an  elliptic  cone  at  several  Mach  number  and  angle 
of  attack  parameter  values.  The  relation  of  the  basic  flows  and  the  linear  disturbances  generated 
upon  them  is  exposed  in  section  4.  The  outlook  with  respect  to  an  inviscid  analysis  of  the  basic 
flows  obtained  is  discussed  in  section  5,  where  the  generalised  Rayleigh  equation  is  presented 
and  discussed.  In  the  light  of  the  interplay  between  basic  states  and  global  linear  disturbances 
alluded  herein  we  present  in  section  6  a  self-contained  discussion  elaborating  on  the  connection 
between  residuals  in  a  direct  numerical  simulation  aiming  at  recovery  of  a  steady-state  solution 
of  the  equations  of  motion  and  global  linear  eigenmodes.  An  efficient  algorithm  for  the  recovery 
of  the  sought  steady-state  using  transient  data  is  presented.  Closing  remarks  are  furnished  in 


section  7. 
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2.  The  Basic  Flow  Computation 

(a)  Governing  equations 

Flow  is  taken  to  be  described  by  the  system  of  the  compressible  continuity,  Navier-Stokes  and 
energy  equations  written  in  conservative  form  as 

dq  a{Fe  +  F4  ,  a{Ge  +  G4  ,  a{He  +  H4 

dt  dx  ^  dy  ^  dz  '  ^  ’ 

with  the  distinction  between  convective  (Fc,  Gc,  Hc)^  and  viscous  (F„G„H,)T  fluxes  being 
made  on  grounds  related  with  the  numerical  discretisation.  The  Cartesian  spatial  coordinates 
are  denoted  by  x,y  and  z  respectively  and  t  denotes  time.  The  vector  of  unknowns  is  q  = 
{p,pu,pv,pw,E)'^,  p  is  the  fluid  density,  {u,v,w)'^  are  the  velocity  components  in  the  directions 
of  X,  y  and  z  and  E  is  the  internal  energy  of  the  fluid.  The  form  of  the  fluxes  is  to  be  found  in 
any  textbook  on  compressible  flow  and  numerical  solutions  thereof  (Hirsch  1988);  the  convective 
fluxes,  of  interest  here  from  the  point  of  view  of  an  inviscid  global  linear  analysis,  are  given  by 


F. 


/  pu 
pu^  +p 
puv 
puw 
V  u{E+p)  ) 


G, 


/  pv 

puv 

pv"^  +p 


(  pw 
puw 
pvw 
pw^  +  p 


(2.2) 


pvw 

\  v{E  +  p)  J  \  w{E+p)  J 

For  the  purposes  of  the  instability  analysis  it  is  convenient  to  recast  the  inviscid  part  of  the 
energy  equation  in  one  of  the  possible  non-conservative  forms 


, du  dv  dw dp  dp  dp.  ^ 

+  fc  >  + ‘“to + 


(2.3) 


One  point  which  deserves  consideration  is  whether  a  perfect  gas  assumption  may  be  made  to 
close  the  system  of  equations;  the  decision  depends  on  the  Mach  number  values  at  which  global 
linear  instability  analysis  of  the  flow  is  to  be  performed  and  affects  the  calculation  of  both  the 
basic  flow  and  the  instability  results.  If  a  perfect  gas  assumption  is  made,  density,  pressure  p, 
temperature  T  and  total  energy  of  the  gas  are  related  by 


'YM'^P  =  pT,  and  p  =  (7  —  1)  ^E  —  (2-4) 

where  7  is  the  ratio  of  specific  heat  coefficients  under  constant  pressure  and  constant  volume, 
7  =  Cp/c^,  and  =  [u^  +v‘^  +  w^)  is  the  kinetic  energy  of  the  fluid  per  unit  mass.  Dimensionless 
numbers  characterising  the  flow  are  the  Mach  number  M,  Reynolds  number  Re  and  Prandtl 
number  Pr,  evaluated  at  appropriate  conditions  and  respectively  defined  as 

M  =  u/og,  Re  =  puLJ p,  and  Pr  =  pcpjk,  (2-5) 

with  u  the  magnitude  of  a  typical  velocity,  Og  the  speed  of  sound,  L  a  typical  length  scale  in 
the  flow,  p  the  dynamic  viscosity  of  the  medium  and  k  the  thermal  conductivity  coefficient.  A 
further  issue  arising  in  hypersonic  flow  is  whether  the  viscosity  may  be  taken  to  be  a  function 
of  temperature  alone  and,  if  so,  what  the  most  appropriate  law  describing  its  dependence  on 
temperature  is.  In  classic  linear  analyses  of  compressible  flow  instability  (Mack  1969,  1984)  a 
power-law  or  Sutherland’s  formula  are  used,  depending  on  whether  the  basic  flow  model  adopted 
is  based  on  a  reduced  or  the  full  system  of  the  equations  of  motion. 
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(b)  Flux- splitting  and  the  AUSMDV  scheme 
A  plethora  of  numerical  schemes  exist  for  the  discretisation  of  (2.1)  most  of  which  use  central 
differencing  for  the  calculation  of  the  viscous  fluxes  (F^,G^,H^)^  and  different  prescriptions 
for  the  calculation  of  the  convective  fluxes  (Fc,  Gc,Hc)^.  A  basic  distinction  exists  between 
shock-capturing  and  shock-fitting  approaches  for  the  discretisation  of  the  convective  fluxes.  We 
discuss  here  briefly  the  scheme  utilised,  namely  the  AUSMDV.  It  should  be  noted  that  a  wealth 
of  modern  accurate  compact  finite-difference  alternative  numerical  methods  exist,  f.e.  the  shock¬ 
capturing  approach  of  Visbal  and  Gaitonde  (1998)  or  the  shock-fitting  algorithm  of  Zhong  (1998, 
1999). 

The  work  of  Liu  and  co-workers  (e.g.  Wada  and  Liu  1994)  presents  the  essentials  of  an  AUS¬ 
MDV  discretisation.  For  the  sake  of  exposition  of  the  algorithm  used  here  we  adhere  with 
presentation  in  one  spatial  dimension.  As  a  first  step,  it  is  recognised  that  the  flux  vector  is 
composed  of  momentum  flux  and  pressure  force  terms,  which  are  treated  independently  as 


Fc  =  u 


+ 


5 


where  E  =  pe.  A  consistent  pressure  splitting  follows,  according  to  which  Fc 


(2.6) 

F+  -|-  F“  with 


F^  =  max(u,  0)  |  pu  j  -|-  j  p“*"  j  ’  ^  ~  min(u,  0)  |  pu  j  -|-  j  p  j  ,  (2-7) 

\pe  J  \  {pu)+  J  \pe  J  \  (pu)-  J 


and 


r  0,  M  <  -1 

p+  =p\  (l  +  M)/2,  \M\  <  1 
[  1,  M  >  1 


r  1,  M  <  -1 

p-=p\  (l-M)/2,  |M|<1  , 

[  0,  M  >  1, 


r  0,  M  <  -1  r  u,  M  <  -1 

(pu)+=p^  {u  +  as)/2,  |M|  <  1  ,  {pu)~  =p\  (u-as)/2,  |M|  <  1 

[  u,  M  >  1  [  0,  M  >  1. 


(2.8) 


The  extension  of  this  idea  to  multiple  spatial  dimensions  is  straightforward.  As  it  stands  this 
scheme  is  only  first-order  accurate  in  space,  while  for  realistic  applications  at  least  second  order 
accuracy  is  necessary  away  from  discontinuities. 


(c)  Numerical  discretisation 

Temporal  discretisation  of  (2.1)  is  performed  using  a  multi-stage  Runge-Kutta  scheme.  For 
the  three-dimensional  problem  at  hand  the  key  additional  issue  arising,  besides  the  choice  of 
a  spatial  discretisation  scheme,  is  that  of  adequate  spatial  discretisation.  This  is  addressed  by 
performing  a  surface  discretisation  of  (1.1)  first,  followed  by  a  consistent  discretisation  of  three- 
dimensional  space  in  which  the  elliptic  cone  is  embedded.  As  a  matter-of-fact  it  has  been  observed 
that  for  the  second-order  scheme  used  best  results  can  be  achieved  when  the  three-dimensional 
spatial  discretisation  is  subdivided  in  two  parts,  first  wrapping  a  structured  grid  around  the 
elliptic  cone,  which  is  generated  conforming  with  the  surface  discretisation  and  then  discretising 
space  between  the  structured  glove  wrapped  around  the  cone  and  the  far-field  boundaries  by  an 
unstructured  grid  generated  by  an  advancing  front  Delauney  algorithm. 

The  grid  used  for  the  a  /  0  calculations  to  be  presented  shortly,  discretising  space  between 
[— ®oo  <  X  <  Xoo]  X  [0  <  y  <  Poo]  X  [—Zoo  ^  z  <  Zoo]  with  Xoo  =  Poo  =  Zoo  =  20,  is  shown  in 
figures  2-4.  The  dimensions  of  the  elliptic  cylinder  are  determined  by  taking  the  aspect  ratio 
to  be  constant  A  =  3  and  the  length  of  the  cone  a  =  1.  The  relevant  parameters  of  Jorgensen 
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(1958)  determine  the  other  two  semiaxes  of  the  cone,  b  =  0.87/3.67  and  c  =  6/3.  In  figure  2 
a  global  view  of  the  discretisation  of  a  half  elliptic  cone  model  surface  and  that  of  the  far- 
field  and  symmetry  boundaries  is  shown,  in  which  the  different  degrees  of  refinement  of  the 
neighbourhood  of  the  object  and  the  far-field  can  also  be  seen.  A  detailed  view  of  the  same  is 
presented  in  figures  3  and  4  viewing  the  half  model  from  upstream  and  downstream,  respectively. 
In  all  three  figures  the  discretisation  of  the  interior  of  the  flowfield  is  not  shown  for  clarity;  in 
the  latter  two  figures  the  structured  grid,  its  fine  distribution  around  the  model  as  well  as  its 
merging  into  the  unstructured  far-field  grid  can  be  seen. 

Thanks  to  the  hybrid-grid  technique  used  the  number  of  total  grid  elements  contained  in  this 
grid  is  confined  to  ~  5.3  x  10®  resulting  in  ~  100  Mbytes  of  memory  required  for  the  storage  of 
the  surface-  and  ~  300  Mbytes  for  that  of  the  field-discretisation,  respectively;  a  much  higher 
number  of  gridpoints  would  have  been  necessary  had  the  structured  grid  in  the  neighbourhood 
of  the  surface  been  extended  to  the  farfield.  The  resolution  and  storage  requirements  increase 
proportionately  to  the  number  of  multigrid  cycles  or  grid  adaptation  levels  in  case  a  multigrid 
or  adaptive-grid  solution  algorithm  is  employed  in  order  to  accelerate  convergence  to  or  refine 
resolution  of  the  steady-state  solution  desired.  We  will  return  to  these  points  when  we  discuss 
basic  flow  results  on  the  elliptic  cone. 


Contract  No.  P61775-99-WE049 


Basic  flows  on  an  elliptic  cone  and  their  global  linear  instability  9 

3.  Basic  Flow  Results 

{a)  Validation  of  the  AUSMDV  scheme 

The  substantial  progress  of  numerical  algorithms  for  external  aerodynamics  over  the  last 
decades  may  be  appreciated  by  referring  to  classic  review  articles  which  periodically  summarise 
the  state-of-the-art  of  the  time.  Such  a  work,  discussing  a  wealth  of  nontrivial  validation  cases, 
is  that  of  Woodward  and  Colella  (1984).  Besides  establishing  the  characteristics  of  the  scheme 
utilised  for  the  recovery  of  the  laminar  steady-state  basic  flows  on  the  elliptic  cone,  we  engage 
in  some  validations  in  order  to  identify  the  limits  of  the  chosen  numerical  approach  on  simple 
configurations;  the  conclusions  drawn  are  then  expected  to  carry  weight  in  the  three-dimensional 
geometry  in  question  on  which  it  is  neither  straightforward  nor  practical  to  perform  a  large 
number  of  numerical  validation  experiments. 

The  solutions  of  two  one-dimensional  Riemann  problems  are  presented  in  figure  5.  The  initial 
conditions  for  the  first  problem  were  proposed  by  Sod  (1978)  and  are  (pl?.^l?Pl)  =  (1?0, 1) 
and  {pr,  Mii,pji)  =  (0.125,0,0.1)  while  those  for  the  second  problem,  due  to  Lax  (1954),  are 
{Pl,  =  (0.445,0.698,^528)  and  {pR,Mji,pji)  =  (0.5,0,0.571).  The  adequate  represen¬ 

tation  of  shocks  and  contact  discontinuities  at  the  resolution  utilised  can  be  seen  in  both  cases. 
A  third  interesting  one-dimensional  problem,  that  of  a  blast  wave,  has  been  discussed  by  Wood¬ 
ward  and  Colella  (1984).  The  periodic  boundary  conditions  of  the  first  two  Riemann  problems 
are  replaced  here  by  reflecting  walls  and  the  calculation  is  initialised  with  the  values  presented  in 
table  1.  Two  strong  shocks  are  generated  and  travel  in  opposite  directions,  are  reflected  by  the 
walls  and  meet  in  the  centre-part  of  the  domain.  The  solution  at  t  =  0.38  is  shown  in  figure  6. 
While  by  visual  comparison  with  the  results  obtained  by  the  piecewise-parabolic  method  of 
Woodward  and  Colella  (1984)  the  solution  obtained  appears  acceptable,  grid  refinement  reveals 
interesting  (but  not  unexpected)  properties  of  the  AUSMDV  method.  Using  current  hardware 
technology  extremely  fine  grids  could  be  used,  the  solution  at  the  aforementioned  time  being 
obtained  within  a  few  minutes  of  CPU  time  on  a  EV6  Alpha  processor.  In  figure  7  the  solution 
at  the  same  time  t  =  0.38  is  presented,  obtained  using  Ax  =  5, 2, 1  and  0.5  x  10“‘*.  In  the  upper 
part  of  this  figure  the  convergence  of  the  solution  at  the  two  highest  resolutions  can  be  seen. 
While  the  location  of  the  shocks  is  adequately  represented  at  all  resolutions,  the  results  at  the 
coarsest  resolution  (which  already  uses  double  the  number  of  points  compared  with  those  pre¬ 
sented  in  figure  6)  is  seen  to  deliver  inaccurate  results  on  at  least  two  counts.  First,  the  details  of 
the  flow  are  smeared  out  in  a  manner  typical  of  numerical  discretisation  schemes  of  low  formal 
order  of  accuracy.  Second,  the  peak  value  of  density  on  the  rightmost  (expansion)  structure  is 
underpredicted,  only  acquiring  its  converged  value  at  the  highest  resolutions.  The  enlargement 
of  the  region  of  the  flow  behind  the  leftmost  shock  reveals  further  details.  Again,  the  coarsest  res¬ 
olution  of  2000  uniformly  distributed  points  smears  the  shock  over  a  rather  wide  distance  while 
as  the  resolution  increases  the  slope  increases  but  pointwise  oscillations  appear  immediately  ad¬ 
jacent  to  the  shock;  their  amplitude  diminishes  with  further  increase  of  resolution.  In  table  2  the 
predicted  wall  values  of  the  primitive  variables  are  presented.  It  may  be  seen  that  even  at  the 
highest  resolution  no  satisfactory  convergence  has  been  achieved,  at  least  from  the  point  of  view 
of  instability  analysis  which  requires  accurate  prediction  of  wall-values  including  their  second 
derivatives.  Of  course,  such  resolutions  are  unattainable  in  the  context  of  the  three-dimensional 
solutions  sought  in  the  framework  of  the  current  investigations  and  a  direct  consequence  of  the 
first-order  accurate  scheme  chosen.  The  conclusions  drawn  from  the  present  one-dimensional 
validation  runs  is  that  the  AUSMDV  method  is  capable  of  predicting  accurately  the  features  of 
unsteady  one-dimensional  flows  encompassing  strong  shocks,  and  is  much  simpler  in  comparison 
with  elaborate  traditional  flux-vector  splitting  schemes,  e.g.  that  due  to  van  Leer.  However,  for 
the  subsequent  instability  analysis  one  should  employ  either  extremely  high  resolutions  (only 
affordable  in  one  or  maximally  two  spatial  dimensions),  or  solution  techniques  based  on  grids 
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other  than  uniform  (f.e.  adaptive  grids).  We  will  return  to  this  subject  shortly,  when  we  present 
the  elliptic  cone  results. 

Before  doing  so,  however,  we  discuss  the  performance  of  the  AUSMDV  scheme  in  two  spatial 
dimensions.  A  plethora  of  nontrivial  two-dimensional  validations  have  been  proposed  in  the 
literature,  of  which  we  have  focussed  our  attention  on  the  two-dimensional  analogon  of  Sod’s 
problem  (in  which  the  initial  conditions  of  the  one-dimensional  problem  are  imposed  in  both 
spatial  directions)  and  the  two  cases  presented  by  Woodward  and  Colella  (1984).  The  result  of 
the  two-dimensional  Sod  problem  may  be  seen  at  time  t  =  2  in  figure  8.  The  interesting  oblique 
pattern  generated  by  the  interaction  of  shocks  and  expansion  waves  and  its  adequate  resolution 
may  be  seen  in  this  result.  However,  this  success  did  not  carry  over  to  the  two-dimensional 
problems  discussed  by  Woodward  and  Colella,  at  least  when  the  first-order  accurate  AUSMDV 
scheme  was  used.  In  both  the  cases  of  flow  over  a  step  and  double  Mach  reflection  difficulties 
with  both  the  accurate  implementation  of  boundary  conditions  and  the  diffusive  nature  of  the 
first-order  accurate  algorithm  have  been  experienced.  While  further  experimentation  has  been 
postponed  until  the  elliptic  cone  flows  have  been  generated,  the  implication  for  the  solutions 
sought  is  clear,  namely  that  a  higher-order  accurate  version  of  the  discretisation  scheme  has  to 
be  utilised  for  the  generation  of  the  elliptic  cone  laminar  solutions;  this  conclusion  is  in  line  with 
that  drawn  earlier  on  the  number  of  discretisation  points  necessary  for  convergence  in  the  one¬ 
dimensional  validation  cases  presented.  With  this  experience  generated,  we  turn  our  attention 
to  the  elliptic  cone  next. 

(6)  Compressible  three-dimensional  laminar  flows  over  an  aspect  ratio  3  elliptic  cone 

In  comparison  with  the  validation  cases  presented  the  calculation  of  flow  over  the  elliptic 
cone  presents  several  additional  challenges,  the  most  obvious  of  which  is  that  related  with  flow 
separation.  A  decision  has  to  be  made  at  the  stage  of  designing  the  computationally  intensive 
solutions  whether  to  embed  the  cone  into  a  calculation  domain  in  which  the  farfield  boundaries 
are  placed  far  away  from  the  three-dimensional  conical  object  or  whether  the  downstream  farfield 
boundary  and  the  surface  of  the  elliptic  cone  meet.  We  have  chosen  the  first  approach,  in 
anticipation  of  inconsistencies  in  the  formulation  of  the  boundary  conditions  at  the  junction  of 
farfield  boundary  and  cone  surface.  While  no  such  problem  is  expected  if  the  downstream  farfield 
boundary  is  placed  several  cone  lengths  downstream  of  the  elliptic  cone,  one  must  decide  on  the 
form  of  the  base  surface;  owing  to  its  simplicity  a  flat  (elliptic)  base  was  chosen.  This  in  itself 
generates  a  new  issue,  namely  unsteady  flow  separation  at  the  base,  which  is  distinct  from  that 
expected  on  account  of  increasing  angle  of  attack  but  will  influence  the  latter  in  the  context  of  a 
three-dimensional  calculation,  at  least  in  the  subsonic  regime.  Since  a  steady  laminar  basic  flow 
is  required  for  the  global  instability  analysis  the  flow  Reynolds  number  must  be  kept  low  enough 
for  such  as  state  to  exist.  To  our  knowledge  the  choice  of  a  Reynolds  number  which  results 
in  steady  separation  on  the  configuration  in  question  is  an  open  question  but  experience  with 
the  classic  Karman  vortex  street  suggests  that  the  Reynolds  number  required  to  achieve  steady 
separation  at  base  may  need  to  be  taken  rather  small.  On  the  other  hand,  given  sufficient  length 
of  the  cone  at  a  moderately  high  Reynolds  number  the  flow  will  become  turbulent  (Kimmel  et 
al.  1997,  1999)  and  a  global  instability  analysis  would  be  obsolete.  With  both  considerations  in 
mind  we  have  taken  a  cone  length  a  =  1  and  Re  =  10^  throughout  the  following  calculations. 
Additionally,  adiabatic  wall  conditions  were  imposed  and  results  were  obtained  in  SI  units. 

In  order  to  avoid  convergence  problems  associated  with  flow  separation  in  the  framework 
of  an  unsteady  calculation,  which  would  be  associated  with  amplification  of  the  global  flow 
eigenmodesf  we  chose  to  perform  steady  basic  flow  calculations.  It  should  be  stressed  at  this 
point  that  the  subsequent  global  instability  analysis  has  the  potential  to  confirm  or  refute  this 


t  see  detailed  discussion  in  §  6 
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choice  in  terms  of  the  damping  or  amplification  of  the  global  eigenmodes  recovered.  The  true 
physical  solution,  at  least  in  its  linear  limit,  can  be  predicted  using  a  two-step  approach;  first 
the  steady-state  equations  of  motion  are  iterated  until  a  steady-state  solution  is  recovered  and 
subsequently  a  global  instability  analysis  of  the  obtained  steady  solution  is  performed.  If  all 
eigenmodes  recovered  by  the  global  analysis  are  damped  the  steady-state  solution  satisfies  the 
unsteady  equations  of  motion  as  well;  if  unstable  eigenmodes  exist  the  computed  steady-state 
solution  will  be  linearly  modified  in  the  manner  predicted  by  the  global  instability  analysis. 

Parameters  which  have  to  be  considered  other  than  the  Reynolds  number  are  the  angles  of 
attack  a  and  bank  fl,  defined  in  the  schematic  representation  of  the  elliptic  cone  geometry 
in  figure  1;  these  determine,  respectively,  flow  separation  from  the  elliptic  cone  surface  and 
symmetry-breaking  of  the  flowfield  around  the  cone.  All  calculations  to  be  presented  have  been 
performed  using  =  0  which  led  to  the  decision  to  impose  symmetry  about  y  =  0.  Whether  this 
decision  is  justified  is  left  to  be  assessed  by  comparisons  with  three-dimensional  solutions  around 
the  entire  elliptic  cone.  Solutions  around  the  entire  cone  are  at  least  twice  as  expensive  as  those 
presented  in  what  follows  and  have  thus  been  left  for  future  extensions/refinements  of  the  present 
work.  Taking  a  =  0  as  well  results  in  the  ability  to  perform  calculations  on  one  quarter  of  the 
elliptic  cone,  imposing  symmetry  about  z  =  0  also.  The  flowfields  resulting  in  this  limit  are  those 
closest  to  fulfilling  the  fundamental  assumption  of  the  global  linear  instability  theory,  according 
to  which  the  flow  dependence  on  two  spatial  directions,  those  defining  the  plane  normal  to  the 
elliptic  cone  surface,  dominate  over  that  along  the  cone  generator.  However,  this  limiting  case 
may  not  be  as  representative  of  realistic  situations  in  which  flow  separation  from  the  elliptic 
cone  surface  on  account  of  a  /  0  and  the  consequent  fulfillment  of  inviscid  instability  analysis 
criteria  is  the  predominant  flow  feature.  In  the  following  calculations  we  have  thus  refrained 
from  imposing  symmetry  about  z  =  0  and  have  considered  two  angles  of  attack,  a  =  10°  and 
20°.  Finally,  the  target  Mach  number  M  of  interest  corresponds  to  hypersonic  flow,  M  =  8.  We 
devised  a  straightforward  but  conservative  strategy  to  achieve  the  latter  target,  commencing 
calculations  at  low  Mach  number  and  angle  of  attack  and  gradually  increasing  the  values  of 
both  parameters  as  a  steady  state  at  a  given  pair  of  (M,  a)  is  approached.  In  this  process  single 
(non-optimal)  grid  calculations  are  performed;  close  to  the  steady  state  a  multigrid  convergence 
acceleration  technique  can  be  employed  and  the  gridpoint  distribution  may  be  optimised  by 
following  the  particular  flow  features  (e.g.  high  gradients).  In  this  respect  a  database  of  steady 
states  at  intermediate  parameter  values  has  been  generated  for  future  reference.  It  should  also 
be  stressed  here  that  at  this  stage  of  the  present  effort  we  are  interested  in  trends  as  parameters 
change  and  not  in  the  machine-accurate  convergence  of  the  solution  at  a  single  set  of  parameters. 
Hence  we  perform  a  large  but  reasonable  number  of  iterations  towards  the  steady  state.  On  the 
hardware  utilised  (EV6  Alpha  processor)  10‘*  iterations  were  performed  at  each  set  of  parameters 
at  a  cost  of  ~  180  CPU  hrs  each.  While  no  statement  can  be  made  regarding  the  convergence  of 
these  results  prior  to  performance  of  simulations  at  the  same  parameters  and  higher  resolutions, 
the  reduction  of  the  residuals  in  the  iterative  solution  by  approximately  three  orders  of  magnitude 
compared  with  the  respective  initial  conditions  warrants  a  qualitative  description  of  the  findings. 

In  figures  9-25  we  present  three-dimensional  flowfields  around  the  aspect  ratio  3  elliptic  cone 
in  the  range  0.5  <  M  <  4.  The  actual  calculation  results  have  been  mirrored  about  y  =  0  in 
order  to  generate  a  three-dimensional  impression.  In  these  figures  the  primitive  variables  in  the 
neighbourhood  of  the  cone  are  shown  as  colour-coded  contour  lines  drawn  on  two  planes,  the 
symmetry  plane  y  =  0  when  available  and  the  plane  meeting  the  cone  at  a:  =  0.7.  We  start  at 
the  subsonic  Mach  number  M  =  0.5  at  which  the  results  of  figures  9-11  and  12-14,  respectively 
pertaining  to  the  angles  of  attack  a  =  10°  and  20°,  are  compared.  Separation  at  the  base  of  the 
cone  is  clearly  visible  in  figures  9  and  11  while  the  contours  of  v  in  figure  10  indicate  vortical 
motion  near  the  surface  of  the  elliptic  cone.  Note  also  that  in  both  these  and  the  results  that 
follow  v{y  =  0)  =  0  on  account  of  the  imposed  symmetry.  However,  it  can  be  seen  in  the  results 
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of  u{z  ^  0)  ^nd  6V6I1  clo&rcr  on  those  of  w(^z  ^  0)  thnt  the  condition  of  symmetry  is  only  satisfied 
in  the  immediate  vicinity  of  the  cone  and  further  iterations  are  necessary  before  the  solution  has 
converged  in  the  far-field  as  well. 

Turning  to  comparisons  with  the  solution  obtained  at  the  larger  angle  of  attack  a  =  20°  one 
notes  that  again,  further  iterations  are  necessary  for  the  condition  of  symmetry  across  y  =  0  to 
be  satisfied.  On  the  other  hand,  comparison  of  the  results  of  figures  10  and  13  gives  a  visual 
impression  of  the  strengthening  of  separation  from  the  cone  surface.  A  better  visualisation  of 
this  result  might  have  been  possible  by  implementing  the  A2  vortex  identification  criterion  of 
Jeong  and  Hussain  (1995);  this  is  currently  in  progress  alongside  detailed  data  analysis.  What 
is  clearly  visible  in  the  contours  presented  is  that  the  separation  zone  at  the  base  of  the  cone  is 
substantially  enlarged  when  the  angle  of  attack  increases.  In  an  incompressible  environment  this 
would  affect  the  flow  on  the  cone  surface  itself;  it  is  expected  that  as  the  Mach  number  increases 
the  effect  of  this  flow  feature,  or  potential  modifications  to  it,  will  not  affect  flow  near  the  tip  of 
the  cone  where  the  global  analysis  is  envisaged  to  be  performed.  Nevertheless,  a  fairing  of  the 
cone  base,  for  instance  by  smoothly  blending  the  cone  itself  with  a  prolate  spheroid/ellipsoid 
body,  might  moderate  the  strength  of  the  separated  area  behind  the  cone;  this  direction  is  also 
currently  being  pursued.  Prom  the  point  of  view  of  the  global  instability  analysis,  it  is  interesting 
to  note  that  the  condition  necessary  for  the  analysis,  namely  that  flow  gradients  along  the  cone 
generators  are  much  smaller  than  those  along  the  wall-normal  direction,  appears  to  be  fulfilled 
in  the  results  at  both  angles  of  attack. 

We  next  keep  the  angle  of  attack  constant  at  a  =  20°  and  increase  the  free-stream  flow  Mach 
number  to  M  =  2;  the  results  for  the  primitive  variables,  obtained  after  a  preset  number  of 
10‘*  iterations,  which  started  using  the  solution  presented  in  figures  12-14  as  initial  condition, 
are  to  be  found  in  figures  15-19.  Prom  a  physical  point  of  view,  the  expected  characteristic  of 
supersonic  flow  Mach  cone  is  to  be  seen  in  these  results.  The  shock  emanates  from  the  cone  tip 
while  expansions  are  to  be  seen  at  the  cone  base  as  a  result  of  the  strong  curvature  jump.  An 
interesting  point  may  be  made  with  respect  to  the  global  analysis  of  this  fiowfield  in  terms  of  its 
potential  self- similarity  with  respect  to  translation  along  the  a:— spatial  direction.  In  figure  20, 
in  which  the  solution  for  w  is  replotted,  one  notices  that  a  transformation  might  exist  within 
a  specific  a:— interval  which  reduces  the  dependence  of  the  basic  flow  from  three  to  two  spatial 
dimensions.  This  issue  is  currently  being  examined;  should  such  a  transformation  be  established 
the  original  objective  of  the  present  effort  to  provide  a  two-dimensional  basic  state  for  the  global 
instability  analysis  would  be  materialised  in  a  manner  consistent  with  the  flow  physics  rather 
by  an  ad-hoc  imposition  of  the  global  instability  analysis  condition  on  the  basic  flow.  Prom  a 
numerical  point  of  view,  the  fact  that  there  is  room  for  improvement  of  resolution  is  also  evident; 
additionally  to  the  boundary  layer  which  must  be  well  resolved  at  all  Mach  numbers  the  shock 
wave  poses  new  challenges  to  the  grid  presented  in  figure  2  in  terms  of  resolution  of  both  the 
neighbourhood  of  the  cone  tip  as  well  as  the  Mach  cone.  One  way  to  improve  resolution  is  by  local 
redistribution  of  the  available  gridpoints  following  some  means  of  grid  adaptation  commensurate 
with  the  gradients  in  the  flow. 

The  need  for  such  a  refinement  of  the  solution  approach  became  stronger  at  the  higher  Mach 
number  of  M  =  4.  While  the  qualitative  features  of  the  flow  at  this  Mach  number  are  analogous 
with  those  at  M  =  2,  the  shock  is  now  stronger  and  thinner.  We  hence  performed  a  large 
number  of  iterations  using  the  solution  at  M  =  2,  a  =  20°  as  initial  condition;  however  the 
rate  of  decay  of  residuals  was  unacceptably  slow  and  we  had  to  resort  to  grid  adaptation, 
increasing  the  number  of  gridpoints  by  30%  compared  with  that  shown  in  figure  2.  Subsequently 
10‘*  more  iterations  were  performed,  producing  the  results  shown  in  figures  21-25.  Qualitatively, 
the  obtained  solutions  appear  to  be  sufficiently  smooth,  though  not  yet  converged  to  machine 
accuracy.  Prom  these  results  it  appears  that  flow  separation  at  the  base  of  the  cone  does  not 
affect  flow  on  the  cone  surface  itself,  the  determining  factors  for  the  latter  being  the  cone  aspect 
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ratio  and  the  angle  of  attack  considered.  However,  as  already  stated,  this  conclusion  is  drawn 
on  the  basis  of  three-dimensional  solutions  of  the  steady  equations  of  motion  and  its  weight  can 
only  be  verified  after  the  global  analysis  has  been  performed.  Our  current  efforts  aim  at  recovery 
of  solutions  at  the  target  Mach  number  M  =  8  using  the  results  at  M  =  4  as  initial  conditions. 
Though  straightforward,  the  achievement  of  this  target  appears  to  be  more  challenging  in  terms 
of  the  computing  effort  required  in  comparison  with  the  solutions  at  lower  Mach  numbers. 
Both  systematic  grid  adaptation  and  multigrid  algorithm  have  been  found  to  be  necessary  in 
preliminary  calculations.  In  parallel,  we  intend  to  pursue  multigrid  solutions  in  order  to  refine 
the  results  at  all  Mach  numbers  shown  herein  in  order  to  achieve  satisfactory  convergence  before 
embarking  upon  data  probing  and  global  instability  analysis. 
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4.  Global  Linear  Instability  Theory 

(a)  Small- amplitude  disturbances  superimposed  upon  a  steady  laminar  flow 

Linear  instability  theory  considers  the  amplification  of  small-amplitude  perturbations  super¬ 
imposed  upon  a  steady  laminar  so-called  basic  flow.  In  its  most  general  form  any  flow  quantity 
is  decomposed  into  an  0(1)  component  denoted  by  and  an  0(e)  perturbation  term  denoted 
by  Qp  according  to 


Q{x,  y,  z,  t)  =  Qb(x,  y,  z)  eQp{x,  y,  z)  exp  Qt  c.c.  (4.1) 

The  vectors  Qfc  =  {p,u,v,w,p)'^  and  Qp  =  {p,u,v,w,p)'^  respectively  denote  the  basic  and 
disturbance  fields,  both  of  which  are  nonperiodic  three-dimensional  functions  of  space.  Taking 
e  <C  1  and  substituting  (4.1)  into  the  governing  system  of  equations  the  0(1)  basic-flow  related 
terms  are  subtracted  out,  themselves  taken  to  satisfy  the  governing  equations  (2.1)  at  this  order. 
Linear  terms  of  0(e)  are  retained  while  higher-order  terms  in  e  are  neglected.  An  eigenmode 
Ansatz  in  time  is  permissible  on  account  of  the  independence  of  the  coefficients  of  the  linearised 
system  of  equations  on  time  and  complex  conjugation  is  considered  in  (4.1)  since  both  and 
the  disturbance  eigenfunctions  Qp  may  in  general  be  complex,  while  both  Q  and  Qj  are  real.  On 
the  other  hand,  the  0(1)  basic  flow  terms,  including  second  spatial  derivatives  thereof,  form  the 
variable  coefficients  of  the  linear  instability  problem;  the  quality  of  the  results  of  the  instability 
problem  critically  depends  on  the  quality  by  which  the  0(1)  terms  are  prescribed.  We  term 
the  analysis  based  on  (4.1)  a  three-dimensional  linear  instability  analysis.  The  study  of  three- 
dimensional  linear  disturbances  satisfying  (4.1)  can  currently  only  be  accomplished  by  a  DNS  for 
Q,  in  which  spatial  derivatives  are  calculated  in  a  decoupled  (usually  parallel)  manner  and  either 
the  linearised  or  the  full  nonlinear  system  of  governing  equations  is  marched  in  time  until  a  regime 
of  exponential  amplification  or  decay  of  perturbations  is  identified  and  results  for  the  temporal 
amplification  rate  of  disturbances  are  obtained  through  dQfdt  in  the  DNS  signal.  A  point 
that  should  be  stressed  here  regards  the  interplay  of  steady  basic  flow  and  three-dimensional 
global  eigenmodes,  namely  that  the  very  existence  of  a  steady  laminar  three-dimensional  basic 
state  Qft  suggests  that  all  three-dimensional  global  linear  disturbances  of  the  flow  are  stable. 
Consequently,  the  performance  of  a  global  linear  analysis  by  numerical  solution  of  the  eigenvalue 
problem,  besides  being  intractable  by  current  hardware  technology,  may  be  of  modest  interest 
from  a  physical  point  of  view;  the  strength  of  a  DNS-based  instability  analyses  then  is  its  ability 
to  address  nonlinear  flow  instability.  Instability  analyses  of  one-dimensional  basic  compressible 
flows  that  utilise  DNS  have  a  somewhat  longer  history  (Erlebacher  and  Hussaini  1990,  Pruett 
and  Zang  1992)  than  linear  analyses  of  two-  and  three-dimensional  basic  states  which  are  slowly 
emerging  (Pruett  and  Chang  1998;  Zhong  1998,  1999). 

Certain  classes  of  basic  flows  may  be  taken  to  be  homogeneous  in  one  spatial  direction  while 
the  other  two  directions  are  resolved.  Specifically,  if 

d/dz  <C  d/dx  and  d/dz  <C  d/dy,  (4.2) 

holds  in  the  basic  flow  problem,  then  the  Ansatz 

Q(a:,  y,  z,  t)  =  Qb{x,  y)  +  eQp{x,  y)  exp  i{A2;  -  Clt}  -\-  c.c.  (4.3) 

can  be  considered  for  linear  disturbances  superimposed  upon  Q^,.  In  this  case  the  linear  in¬ 
stability  analysis  may  proceed  either  by  means  of  DNS  in  a  manner  conceptually  analogous 
with  that  used  for  the  instability  study  based  on  (4.1)  or  by  recasting  the  linearised  system 
of  equations  as  a  partial-differential-equation  based  eigenvalue  problem.  In  a  temporal  frame¬ 
work  the  objective  becomes  the  determination  of  and  the  associated  amplitude  eigenfunctions 
Qp  =  (p,  fi,  V,  w,p)'^  given  a  periodicity  length  in  the  homogeneous  z  spatial  direction  which  is 
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associated  with  the  real  wavenumber  parameter  A  through  A  =  27^  Here  Qp  =  (p,  u,  v,  w,p)'^ 
are  again  three-dimensional  functions  of  space  but  in  contrast  to  those  in  (4.1)  they  are  periodic 
in  z.  Prerequisite  for  the  analysis  in  either  a  DNS  or  an  eigenvalue  problem  framework  is  the 
existence  of  a  steady  (or  time-periodic)  laminar  two-dimensional  state  Qb{x,y).  The  results  of 
the  two  approaches  are  identical  as  far  as  the  leading  instability  eigenmode  is  concerned,  while 
the  eigenvalue  problem  has  the  additional  merit  over  DNS  of  being  able  to  provide  details  of  the 
flow  eigenspectrum  which  are  difficult  either  to  obtain  or  to  identify  in  the  DNS.  We  term  the 
analysis  based  on  Ansatz  (4.3)  a  two-dimensional  linear  instability  approach,  while  on  account 
of  the  inhomogeneity  of  more  than  one  spatial  directions  in  both  the  basic  flow  problem  and 
that  of  the  corresponding  linear  theory,  the  instability  theories  based  on  either  (4.1)  or  (4.3)  are 
termed  global  linear  analyses. 

By  contrast  to  the  global  analyses,  consideration  of  two  spatial  directions  as  homogeneous  and 
their  consequent  resolution  by  Fourier  series  leads  to  the  classic  local  linear  instability  theory,  in 
which  the  so-called  parallel-flow  approximation  is  made  for  the  basic  flow,  such  that 

dQb/dx^d/dy,  dQb/dz^d/dy.  (4.4) 

The  corresponding  decomposition  becomes 

Q(a:,  y,  z,  t)  =  Qb(y)  +  eQp(y)  exp  i{KX  -h  Az  -  fit}  -h  c.c.  (4.5) 

We  term  the  local  linear  analysis  based  on  (4.5)  an  one- dimensional  analysis.  The  analogon 
of  Ansatz  (4.5)  in  a  circular  cone  geometry  is 

Q(r,  (9,  z,  t)  =  Qb{r)  +  eQp(r)  exp  i{n6  -\-  \z  —  fit}  c.c.  (4-6) 

where  r  denotes  the  radial,  6  the  azimuthal  and  z  the  downstream  direction  and  n  is  taken 
to  be  an  integer.  While  a  DNS  is  also  possible  and  indeed  common  practice  as  far  as  nonlinear 
instability  analysis  of  compressible  flows  is  concerned  (e.g.  Pruett  and  Zang  1992)  an  one¬ 
dimensional  linear  instability  analysis  is  efficiently  performed  by  means  of  numerical  solution  of 
the  ordinary-differential-equation  based  generalised  eigenvalue  problem  which  results  when  (4.5) 
is  substituted  into  the  governing  equations,  the  resulting  system  is  linearised  and  terms  of  0(e) 
only  are  retained  and  solved  for. 

A  local  linear  instability  analysis  is  the  first  to  be  considered  from  a  physical  point  of  view 
and  the  most  straightforward  to  perform  as  far  as  numerical  solutions  of  the  eigenvalue  problem 
are  concerned.  Mack  (1969,  1984)  presented  excellent  summaries  of  the  early  work  which  led 
to  the  discovery  of  instability  modes  particular  to  compressible  flow  while  Malik  (1991)  has 
discussed  alternative  numerical  methods  for  the  accurate  and  efficient  solution  of  the  ODE- 
based  linear  eigenvalue  problem.  Owing  to  the  numerical  challenges  presented  by  the  global 
linear  analyses,  the  vast  majority  of  linear  instability  work  performed  to-date  is  confined  within 
the  framework  of  the  local  approach  based  on  (4.5).  The  one-dimensional  basic  flow  Qb(y)  is 
obtained  using  approaches  of  different  degrees  of  sophistication,  based  on  numerical  solutions  of 
the  boundary-layer,  the  Thin-Layer-  or  the  Parabolised-Navier- Stokes  equations  (e.g.  Kimmel  et 
al.  1999).  Within  the  framework  of  an  inviscid  one-dimensional  linear  instability  analysis  either  a 
consistent  approach  is  taken,  typically  in  open  systems,  in  which  an  inviscid  steady  laminar  basic 
flow  forms  the  variable  coefficients  of  the  linear  eigenvalue  problem  or  a  non-rational  approach 
is  followed,  typically  in  wall-bounded  flows,  in  which  viscosity  is  retained  in  the  basic  but  not 
in  the  disturbance  equations.  On  the  other  hand,  within  the  framework  of  the  two-dimensional 
linear  instability  analysis  which  forms  the  target  of  the  present  investigations,  a  steady  laminar 
two-dimensional  flow,  i.e.  a  basic  flow  consisting  of  all  three  velocity  components,  density  and 
pressure,  all  of  which  are  functions  of  two  spatial  coordinates,  must  be  provided. 
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( b )  Comments  on  the  quality  of  the  steady  laminar  basic  flows 

As  far  as  the  basic  flows  Qj  necessary  for  the  analyses  are  concerned,  short  of  resorting  to 
DNS,  these  are  obtained  by  means  of  approaches  consistent  with  the  subsequent  linear  instability 
analysis  performed.  Boundary  layer  (Mack  1969,  Duck  1990,  Karabis  et  al.  1999)  and  Parabolised 
Navier-Stokes  (PNS)  solutions  (Kimmel  et  al.  1997,  1999)  have  been  obtained,  as  well  as  solutions 
delivered  by  Navier-Stokes  solvers  appropriate  for  aerodynamic  calculations  (Dietz  and  Hein 
1999).  In  all  cases  crucial  for  the  success  of  the  subsequent  instability  analysis  is  the  accuracy  by 
which  the  flowfield  is  described.  Specifically,  second  derivatives  of  flow  quantities  in  the  radial 
direction  must  be  accurately  provided  in  the  basic  flow,  a  result  straightforward  to  achieve  in  the 
context  of  boundary-layer  but  nontrivial  in  the  context  of  PNS  or  Navier-Stokes  computations. 
Indeed,  in  the  elliptic  cone  problem  of  interest  here,  the  basic  flow  issue  is  far  from  having  been 
resolved  in  a  satisfactory  manner,  with  not  only  different  approaches  (e.g.  boundary-layer  vs. 
Navier-Stokes)  but  also  different  calculations  within  the  same  approach  delivering  results  for  the 
same  configuration  which  are  in  substantial  disagreement  with  each  other  (Kimmel  et  al.  1997). 
The  implications  for  any  type  of  subsequent  linear  instability  analysis  of  discrepancies  in  the 
basic  flow  are  evident. 

Furthermore,  regarding  a  global  instability  analysis  an  additional  constraint  is  posed  by  the 
number  of  discrete  points  on  which  basic  flow  information  can  be  provided  so  that  the  instabil¬ 
ity  analysis  may  be  performed  within  computationally  affordable  limits.  In  other  words,  what 
is  necessary  is  a  basic  field  Qj  of  highest  quality  on  as  low  a  number  of  points  as  possible.  This 
requirement  is  typically  fulfilled  by  employing  numerical  methods  of  high  formal  accuracy,  typ¬ 
ically  spectral  (Canute  et  al.  1987)  or  compact  finite-difference  schemes  (Visbal  and  Gaitonde 
1998;  Gaitonde  and  Visbal  1999).  Experience  with  both  instability  analyses  and  DNS  comparing 
spectral  and  high-order  compact  finite-difference  methods  has  shown  that  either  discretisation 
delivers  satisfactory  results,  although  the  compact  finite-difference  methods  methods  typically 
require  three-times  the  number  of  discretisation  points  per  spatial  direction  in  order  to  match 
the  accuracy  of  the  spectral  methods  (Theofilis  19986).  In  the  context  of  the  coupled  spatial 
discretisation  of  two  spatial  directions  a  spectral  method  may  thus  be  preferable. 

However,  in  compressible  flow  the  issue  of  existence  and  location  of  shocks  must  be  considered 
in  conjunction  with  linear  instability  analysis.  Instability  analyses  based  on  the  decomposition 
(4.5)  have  to-date  typically  been  performed  on  the  assumption  of  exclusion  of  shocks  from 
the  computation  domain.  This  assumption  must  be  viewed  critically  in  the  context  of  one¬ 
dimensional  compressible  linear  analyses,  where  two  types  of  instabilities  are  known  to  exist,  the 
compressible  analoga  of  the  Tollmien-Schlichting  instabilities  of  incompressible  flow  (Tollmien 
1929)  and  so-called  Mack- modes  particular  to  compressible  flow  (Mack  1969).  Both  types  of 
eigenmodes  decay  rapidly  outside  the  boundary  layer  but  at  a  different  rate,  with  the  Mack 
modes  penetrating  deeper  into  the  free-stream  than  the  TS-like  instabilities.  In  the  case  of  conical 
geometries  inviscid  analysis  has  revealed  the  existence  of  instabilities  the  eigenfunctions  of  which 
are  of  oscillatory  nature  with  their  rate  of  decay  being  a  nonlinear  function  of  the  respective 
eigenvalue  (Duck  1990,  Shaw  and  Duck  1992).  The  consistency  of  the  results  of  such  analyses 
has  to  be  verified  a  posteriori  by  comparison  with  the  shock  location  at  specific  flow  conditions. 
For  both  local  and  global  linear  analyses  the  condition  must  be  fulfilled  that  the  shock  position 
is  well  beyond  the  boundary-layer  edge,  such  that  either  homogeneity  of  the  perturbations  or 
asymptotic  boundary  conditions  based  on  the  existence  of  realisable  free-stream  conditions  may 
be  imposed.  These  considerations  have  led  to  the  identification  of  the  two  alternatives  for  the 
calculation  of  laminar  steady-state  basic  flows  in  the  present  work,  a  shock-capturing  scheme 
based  on  one  variant  of  the  AUSMDV  method  (Wada  and  Liu  1994)  and  a  shock-fitting  approach 
based  on  DNS  (Zhong  1998).  On  account  of  satisfactory  experience  gathered  with  the  former 
scheme  on  a  related  problem  (Hein  et  al.  1999)  the  AUSMDV  scheme  was  chosen. 
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5.  The  Generalised  Rayleigh  Equation 


(a)  Two-dimensional  inviscid  global  linear  analysis 

Returning  to  the  global  instability  problem,  substitution  of  (4.1)  into  (2.1)  leads  to  the  lin¬ 
earised  three-dimensional  compressible  continuity,  Navier-Stokes  and  energy  equations.  In  super¬ 
sonic  flow  Mack  (1969)  has  shown  that  the  essential  physical  one-dimensional  linear  instability 
mechanisms  may  be  well-understood  from  an  inviscid  analysis  point  of  view,  viscosity  moder¬ 
ately  modifying  the  quantitative  results  but  not  their  qualitative  trends.  Here  we  discuss  the 
inviscid  linearised  system  first,  as  well  as  possible  simplifications  which  make  the  global  linear 
instability  problem  tractable  from  a  numerical  point  of  view.  The  assumption  of  a  steady  basic 
flow 


dp  du  dv  dw  _ 

dt  dt  dt  dt 


(5.1) 


leads  to  the  three-dimensional  inviscid  linearised  disturbance  equations,  in  symbolic  form 
written  as 

dQp 


•^SdQp  —  ^Zd- 


dt 


(5.2) 


As  has  been  mentioned,  the  discretisation  of  the  linear  operators  Azd  and  B^d-,  in  which 
the  three-dimensional  basic  flow  is  used  as  the  variable  coefficients,  by  means  of  an  eigenvalue 
problem  in  which  all  three  spatial  directions  are  treated  in  a  coupled  manner  is  impractical  with 
current  hardware  technology  and  uninteresting  from  a  physical  point  of  view,  since  the  existence 
of  a  converged  three-dimensional  steady-state  solution  is  synonymous  with  stability  of  all  three- 
dimensional  global  flow  eigenmodes.  In  order  to  proceed  with  the  analysis  the  simplification  is 
made  that  the  basic  flow  is  independent  of  one  spatial  direction,  say  z,  f  i.e. 

dp  du  dv  dw  dp  _ 

This  permits  the  introduction  in  (5.2)  of  eigenmodes  for  the  disturbance  quantities  in  z  and 
t,  such  that 


Qp{x,  y,  z,  t)  =  Q{x,  y)  exp  i(A2;  -  fit) 


and  results  in  a  linearised  system  of  the  form 


«4.2dQ  —  fiH2<iQ- 


(5.4) 


(5.5) 


where  Q  =  (p,  u,  f),  w,p)'^  are  the  two-dimensional  amplitude  functions  of  the  three-dimensional 
perturbations  Qp{x,y,  z,  t).  The  entries  of  .42^  and  B2d  may  be  found  in  Appendix  A.  As  it 
stands,  the  linearised  system  (5.5)  offers  no  numerical  simplification  over  the  incompressible 
primitive- variables  formulation  discussed  by  Theofilis  (1998c);  indeed,  the  existence  of  an  ad¬ 
ditional  equation  further  aggravates  the  already  large  memory  and  runtime  requirements  for 
numerical  solution  of  the  partial-differential-equation  based  eigenvalue  problem.  If  a  complex 
shift-and-invert  is  performed  in  order  for  one,  instead  of  two  matrices  appearing  in  (5.5)  to  be 
stored,  in  excess  of  0(1.5Gbytes)  of  storage  for  the  resulting  matrix  are  necessary;  an  interesting 
window  of  its  eigenvalues  may  be  calculated  by  Krylov  subspace  iteration  methods  at  a  cost  of 
approximately  one  CPU  hour  of  supercomputing  time  at  a  single  flow  condition.  As  a  matter  of 
fact,  the  solution  of  the  compressible  viscous  global  linear  instability  eigenvalue  problem  requires 
the  same  amount  of  computing  effort  as  the  inviscid  problem  (5.5)  so  that  solution  of  the  viscous 
problem  may  be  preferable  from  the  point  of  view  of  consistency  of  the  basic  flow  and  the  global 
instability  analysis  approaches. 


t  Attention  is  drawn  to  the  fact  that  the  homogeneous  direction  is  denoted  by  z  in  this  section,  and  does  not  refer  /  is 
not  to  be  confused  with  the  coordinate  system  used  on  the  elliptic  cone 
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(b)  Further  simplifications 

Several  simplifications  of  the  system  of  inviscid  linear  disturbance  equations  are  offered  by  a 
basic  flow  which  possesses  a  single  velocity  component,  that  in  the  direction  of  the  wavenumber 
vector.  Physically,  this  situation  is  realisable  by  considering  zero-angle-of-attack  flow  over  the 
elliptic  cone.  The  system  (5.5)  then  becomes 


iXwp  +  {pDx  +  Pxl)u  +  {pDy  +  Pyl)v  +  iXfiw  =  iQp  (5.6) 

iXpwu  +Px  =  iflpu  (5-7) 

iXpwv  +Py  =  iflpD  (5-8) 

pWxU  +  pWyV  +  iXpww  +  iAp  =  iQpw  (5-9) 

{jpDx  +  Pxl)u  +  {"ypDy  +  Pyl)v  +  iXjpw  +  iXwp  =  iClp.  (5.10) 


A  first  simplification  offered  by  (6.9-5.10)  in  comparison  with  (5.5)  is  that  the  former  may  be 
written  as  a  real  eigenvalue  problem,  requiring  half  the  storage  and  computing  time  to  solve  in 
comparison  with  the  latter  problem.  This  fact  may  be  seen  clearly  by  eliminating  p,  u,  v  and  w 
from  (6.9-5.10),  which  results  in 


Cp  + 


“^Xwx 
[Aid  —  fl] 


Px  + 


2XWy 

[Au)  —  fl] 


Py  + 


( p[Xw  — 

V  7P  / 


p  =  0  (5.11) 


where  £  =  Dxx  +Dyy  —  X^,Dx  =  djdx  and  Dy  =  djdy.  This  is  the  compressible  counterpart 
of  the  Rayleigh  equation  within  the  framework  of  two-dimensional  linear  instability,  in  which 
the  (Mach-number  related)  terms  arising  from  compressibility  are  grouped  in  brackets.  The 
incompressible  analogon  of  (5.11)  has  been  discussed  by  Hall  and  Horseman  (1991)  in  the  context 
of  inviscid  secondary  instability  of  Gdrtler  vortices. 

Prom  a  numerical  point  of  view  it  may  be  seen  that  the  major  difference  between  the  incom¬ 
pressible  and  the  compressible  two-dimensional  Rayleigh  equations  is  that  the  former  constitutes 
a  linear  while  the  latter  is  a  cubic  eigenvalue  problem  in  (in  temporal)  or  in  A  (in  spatial) 
global  linear  theory.  The  nonlinearity  in  the  eigenvalue  might  lead  to  using  an  iterative  method 
for  the  solution  of  (5.11).  In  classic  one-dimensional  inviscid  linear  instability  analyses  (Mack 
1969)  one  deforms  the  path  of  integration  into  the  complex  plane  with  a  uniquely  determined 
contour  indentation;  however,  this  uniqueness  is  lost  in  the  context  of  inviscid  global  linear 
calculations  and  the  issue  of  the  correct  integration  path  in  analyses  in  which  two  complex  spa¬ 
tial  coordinates  have  to  be  taken  simultaneously  into  account  is  far  from  having  been  resolved 
(Duck,  personal  communication).  A  direct  solution  algorithm,  on  the  other  hand,  requires  use  of 
a  companion  matrix  approach  (Bridges  and  Morris  1984;  Theofilis  1995)  which  results  in  a  dis¬ 
crete  problem  whose  size  equals  that  of  numerical  solution  of  three  equations.  This  is  the  second 
simplification  offered  by  the  two-dimensional  Rayleigh  equation  in  comparison  with  (6.9-5.10), 
namely  the  savings  of  two  out  of  five  equations,  which  results  in  a  compressible  two-dimensional 
global  eigenvalue  problem  which  is  actually  less  expensive  than  viscous  incompressible  global 
linear  instability  in  primitive  variables  (Theofilis  1998c). 

The  third  simplification  that  the  existence  of  a  single  basic  velocity  component  permits  is 
the  possibility  to  impose  symmetries  in  the  seeked  disturbance  solutions.  Symmetries  have  been 
imposed  in  order  to  tackle  the  incompressible  two-dimensional  linear  instability  problems  in 
a  rectangular  duct  (Tatsumi  and  Yoshimura  1990)  and  in  the  infinite  swept  attachment-line 
boundary  layer  (Lin  and  Malik  1996).  In  the  case  of  hypersonic  flow  over  an  elliptic  cone  the 
instability  problem  could  be  addressed  in  the  context  of  numerical  solutions  in  one  quadrant 
alone,  resulting  in  a  matrix  of  size  sixteen  times  smaller  than  that  quoted  earlier  as  necessary 
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for  solution  of  the  incompressible  global  linear  instability  problem  in  primitive  variables.  Such 
savings  offers  much  more  flexibility  regarding  good  resolution  of  the  global  eigenfunctions  in 
comparison  with  the  full  problem  based  on  (5.5)  and  greatly  facilitates  parametric  studies, 
albeit  only  in  one  of  the  physically  interesting  cases,  that  of  zero-angle  of  attack  flow.  In  order 
to  address  flow  oncoming  obliquely  to  the  axis  of  symmetry  of  the  elliptic  cone  the  full  problem 
(5.5)  must  be  addressed  and  this  is  done  best  in  the  context  of  a  consistent  viscous  calculation. 
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6.  On  numerical  residuals  and  physical  instabilities  in  steady-state  fluid  flow 

calculations 

The  preceding  discussion  repeatedly  refers  to  the  interplay  between  the  existence  of  a  steady- 
state  solution  of  the  equations  of  motion  and  amplification  of  global  eigenmodes  of  this  solution. 
Our  concern  here  is  with  the  identification  of  the  origin  and  classification  of  the  different  qualita¬ 
tive  forms  which  the  numerically  obtained  transient  solution  assumes,  when  a  DNS  is  performed 
in  order  to  recover  a  steady-state  solution  of  the  system  of  equations  governing  fluid  flow  motion. 
For  reasons  of  feasibility  of  the  related  partial- derivative  eigenvalue  problem  we  restrict  the  dis¬ 
cussion  to  the  recovery  of  two-dimensional  steady  basic  states.  Within  the  unifying  framework  of 
the  extension  of  the  linear  instability  theory  of  Tollmien  (1929)  to  nonparallel  two-dimensional 
steady  basic  states,  residuals  encountered  in  the  simulation  as  the  latter  approaches  convergence 
are  either  identified  as  or  associated  with  the  least  damped  of  the  two-dimensional  global  linear 
eigenmodes  of  the  steady-state  flow.  The  inability  to  converge  to  a  steady-state  is  shown  to  be 
linked  with  the  global  linear  flow  eigenmodes  approaching  a  neutrally  stable  state  and  interact¬ 
ing  nonlinearly.  With  the  origin  of  the  residuals  established,  an  algorithm  is  presented  which 
permits  recovery  of  the  converged  steady-state  solution  from  transient  data  at  substantially  less 
computing  effort  compared  with  that  necessary  for  the  integration  of  the  system  of  equations 
until  convergence.  Nonparallel  linear  instability  theory  of  the  three-dimensional  eigenmodes  of 
the  converged  two-dimensional  steady-state  may  also  be  used  to  quantify  the  differences  between 
the  results  of  two-  and  three-dimensional  DNS.  Rather  than  remaining  within  the  framework  of 
the  elliptic  cone  we  use  the  well- documented  two-dimensional  incompressible  lid-driven  cavity 
flow  as  a  demonstrator  of  these  ideas. 

(a)  The  problemacy  with  residuals 

A  steady-state  solution  Qj  of  the  two-dimensional  incompressible  continuity  and  Navier-Stokes 
equations  which  describe  flow  in  a  prescribed  two-dimensional  domain  S  bounded  by  dT,  is  sought 
numerically.  A  plethora  of  numerical  approaches  for  the  accurate  and  efficient  integration  of 
either  the  steady  or  the  unsteady  equations  of  motion  exists  (e.g.  Kim  and  Moin  1985)  so  that 
this  problem  may  be  considered  solved  in  principle.  However,  in  performing  a  time-accurate 
integration  of  the  equations  of  motion  one  observes  that,  depending  on  the  values  of  parameters 
such  as  the  flow  Reynolds  number,  in  the  limit  of  large  time  either  a  steady-state  solution  is 
obtained  (e.g.  Briley  1971)  or  unsteady,  sometimes  periodic,  motion  sets  in  (e.g.  Pauley  et  al. 
1990;  Goodrich  et  al.  al.  1990).  The  first  question  arising  is  what  type  of  physical  information  is 
not  considered  by  solving  the  steady  as  opposed  to  the  unsteady  equations  of  motion  and  what 
is  the  physical  interpretation  of  the  critical  conditions  beyond  which  the  steady  and  unsteady 
formulations  deliver  different  results.  Both  physical  and  numerical  experience  suggest  that  at  low 
Reynolds  numbers  the  two  formulations  may  be  used  interchangeably.  For  example,  essentially 
identical  results  with  those  of  Ghia  et  al.  (1982)  and  Schreiber  and  Keller  (1983)  have  been 
obtained  by  a  multitude  of  subsequent  investigators  who  used  the  time-dependent  equations  of 
motion  to  describe  flow  in  the  square  lid-driven  cavity  at  Reynolds  numbers  up  to  Re  =  10‘*.  On 
the  other  hand,  the  question  of  existence  of  a  steady-state  solution  delivered  by  the  unsteady 
version  of  the  equations  of  motion  at  Re  =  10‘*  has  been  recently  re-opened  (E  and  Liu  1996) 
while  it  is  well  known  that  Hopf  bifurcations  exist  in  both  the  aspect-ratio  two  singular  lid- 
driven  cavity  at  Re  <  5000  (Goodrich  et  al.  1990)  and  its  regularised  square  counterpart  at 
Re  «  10‘*  (Shen  1991).  Consensus  exists  that  at  high  Reynolds  numbers  the  unsteady  formulation 
is  capable  of  delivering  physics  inaccessible  to  the  steady  version  of  the  equations  of  motion; 
however,  the  origin  of  the  differences  between  the  results  of  the  two  formulations  is  presently 
not  understood  in  a  satisfactory  manner.  This  is  an  alternative  way  of  posing  the  first  question, 
namely  what  are  the  unsteady  effects  that  manifest  themselves  at  high  Reynolds  numbers? 

The  next  question  arises  from  the  very  concept  of  two-dimensionality.  The  results  of  numerical 
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solutions  of  the  three-dimensional  analoga  of  the  incompressible  continuity  and  Navier-Stokes 
equations  are  in  most  cases  in  substantial  qualitative  and  quantitative  disagreement  with  their 
two-dimensional  counterparts  (Burggraf  1966)  relegating  two-dimensional  DNS  to  the  realm 
of  academic  interest.  Within  the  scope  of  two-dimensional  solutions  being  of  interest  three- 
dimensionality  of  physical  space  could  be  addressed  by  considering  the  flow  to  be  independent 
of  the  third  spatial  direction.  Homogeneity  in  this  third  direction  could,  in  turn,  be  discussed  in 
the  context  of  a  three-dimensional  simulation,  nonperiodic  in  the  same  two  spatial  directions  as 
the  two-dimensional  one  and  periodic  in  the  third.  Advances  in  both  algorithms  and  hardware 
and,  not  least,  a  considerable  amount  of  knowledge  on  the  differences  between  two-  and  three- 
dimensional  numerical  simulation  results  lead  one  to  employ  a  DNS  algorithm  for  flow  with  two 
nonperiodic  and  one  periodic  spatial  direction  (e.g.  Spalart  1988)  in  the  founded  expectation  that 
a  three-dimensional  so-called  ’spatial’  DNS  is  the  only  means  capable  of  capturing  all  physical 
phenomena  at  a  certain  Reynolds  number.  The  second  question  which  may  be  posed  at  this  point 
regards  the  origin  of  the  differences  between  the  results  of  such  two-  and  three-dimensional  direct 
numerical  simulations.  Associated,  one  may  ask  whether  there  exists  an  alternative  means  to 
spatial  DNS  for  the  description  of  the  origins  of  the  three-dimensional  phenomena  encountered. 

The  objective  of  the  present  effort  is  to  put  both  questions  within  the  unified  framework  of 
nonparallel  linear  instability  of  the  steady  state  Qj,.  With  the  aid  of  the  well-studied  lid-driven 
cavity  flow  example  we  demonstrate  the  intimate  link  between  numerical  residuals  in  steady- 
state  fluid  flow  calculations  and  linear  two-dimensional  eigenmodes  of  the  converged  steady  state 
Qj.  In  the  next  section  6  (b )  we  present  theoretical  arguments,  first  analysing  the  behaviour  of 
numerical  residuals  near  convergence  towards  the  steady-state  solution  from  a  numerical  point 
of  view.  Subsequently  we  discuss  solutions  of  the  partial  derivative  eigenvalue  problem  governing 
linear  instability  of  nonparallel  two-dimensional  steady-state  flows,  which  shed  light  on  residuals 
from  a  physical  viewpoint.  With  the  origin  of  residuals  established  from  a  physical  point  of  view 
we  construct  and  present  an  algorithm  which  permits  recovery  of  the  converged  steady-state 
solution  from  transient  results  of  the  time-marching  procedure,  the  latter  taken  well  before 
convergence.  In  section  6  (c )  we  present  results  of  application  of  the  algorithm.  First,  the  link 
between  the  results  of  nonparallel  linear  instability  theory  and  different  types  of  behaviour  of 
numerical  residuals  in  the  DNS  is  demonstrated  in  this  section  and  the  aforementioned  questions 
are  answered.  Subsequently,  examples  of  recovery  of  the  converged  steady-state  from  transient 
data  and  assessment  of  the  substantial  savings  in  the  computing  effort  materialised  by  use 
of  the  proposed  algorithm  are  presented  in  this  section.  Closing  remarks  on  the  far-reaching 
implications  of  the  present  findings  are  made  in  section  6  (d  )  where  suggestions  for  the  extension 
of  the  analysis  presented  to  compressible  flow  and  flow  with  three  nonperiodic  spatial  directions, 
such  as  that  on  the  elliptic  cone,  are  made. 

(6)  Theory 

(i)  On  residuals  and  the  phenomenology  oe  their  behaviour 

While  in  an  computation  based  on  the  steady  system  of  equations  governing  fluid  flow  motion, 
such  as  that  used  for  the  basic  flows  on  the  elliptic  cone  presented  in  the  previous  section, 
residuals  are  viewed  as  departure  from  the  steady  state  which  have  to  be  eliminated  in  an 
efficient  manner  by  a  specific  solution  algorithm  (e.g.  multigrid),  in  a  time-accurate  integration 
one  may  view  transients  as  solutions  of  the  equations  of  motion  at  every  time-step  and  attempt 
to  attach  physical  significance  to  characteristic  patterns  of  their  behaviour.  By  contrast  to  the 
previous  sections  here  we  concentrate  on  a  time-accurate  integration  of  the  unsteady  equations 
of  motion  and  monitor  the  behaviour  of  residuals,  defined  as  the  difference  between  the  transient 
solution  and  the  converged  steady  state,  in  flow  regimes  where  the  latter  exists.  Physical  space 
is  three-dimensional;  without  loss  of  generality  we  may  take  the  Cartesian  coordinates  x  and  y 
to  be  defined  on  S  while  z  denotes  the  third  spatial  coordinate  in  the  direction  of  S.  Along  the 
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first  two  coordinates  the  velocity  vector  has  components  u  and  v,  while  pressure  is  denoted  by 
p.  The  equations  of  motion  are  marched  in  time  t  until  q  =  {u,v,p)'^,  the  transient  solution, 
converges  to  Q^.  Assuming  that  the  latter  exists  and  keeping  the  domain  S  unchanged,  the 
following  qualitative  observations  are  made. 

First,  at  any  Reynolds  number  Re  at  which  Qj  exists,  close  to  convergence  the  residuals 
decay  exponentially  in  amplitude.  Second,  refinement  of  the  discretisation  of  the  domain  S  at 
constant  Re  results  in  convergence  of  the  rate  at  which  the  residuals  decay.  Third,  the  (converged) 
rate  of  decay  of  residuals  is  a  function  of  the  flow  Reynolds  number;  as  Re  increases  residuals 
decay  slower  and  the  associated  time  of  integration  of  the  equations  of  motion  until  convergence 
increases.  Fourth,  on  occasion,  the  residuals  decay  at  a  specific  constant  rate  for  a  number  of 
decades  before  this  rate  of  decay  changes  to  a  different  constant  value  at  which  residuals  further 
decay  until  convergence.  Fifth,  systematically  increasing  Re,  instead  of  monotonic  convergence 
of  residuals  an  oscillatory  behaviour  of  q  in  the  neighbourhood  of  is  observed.  Ultimately,  a 
value  of  Reynolds  number  is  reached  past  which  no  Qj  exists.  At  first  sight  the  existence  of  a 
physical  mechanism  which  unifies  such  diverse  patterns  of  behaviour  of  the  numerical  solution 
seems  unlikely. 

(ii)  A  NUMERICAL  POINT  OF  VIEW  ON  THE  BEHAVIOUR  OF  RESIDUALS  NEAR  CONVERGENCE 

However,  it  is  straightforward  to  provide  an  explanation  of  the  first  observation  on  the  be¬ 
haviour  of  residuals,  which  also  provides  a  handle  to  the  link  between  numerical  residuals  and 
physical  flow  instabilities.  We  assume  that  the  solution  q  is  close  to  converging  to  the  sought 
two-dimensional  field  Qj  =  {u,v,p)'^  such  that  it  may  be  decomposed  into  the  latter  and  small 
two-dimensional  residuals  q2D  =  {u2i),V2i3,P2i))'^  superimposed  upon  it,  according  to 

q(a:, y, t)  =  Qb{x, y)+e  ci2T>{x, y, t),  (6.1) 

with  e  <C  1.  We  next  substitute  the  decomposition  (6.1)  into  the  continuity  and  Navier-Stokes 
equations  and  assume  that  the  steady-state  solution  satisfies  the  equations  of  motion  at  0(1), 
such  that  it  may  be  subtracted  out  of  the  resulting  system  at  this  order.  Subsequently,  based 
on  the  smallness  of  the  amplitude  of  the  residuals,  we  linearise  about  and  rearrange  the 
system  at  0(e)  such  that  the  vector  of  residuals  represents  the  unknowns;  terms  of  O(e^)  are 
neglected.  Since  the  coefficients  of  the  resulting  linear  system  of  equations  for  the  determination 
of  q2D  at  0(e)  are  independent  of  time  t  we  may  introduce  an  eigenmode  decomposition  in  this 
coordinate,  according  to 

h2T>{x,y,t)  =  q2D(a:,y)  (6.2) 

with  q2D  =  {u2T),V2T),P2T))'^ ■  The  physical  significance  of  the  parameter  a  will  be  discussed 
shortly;  from  a  numerical  point  of  view  it  represents  the  rate  at  which  the  residuals  q2D  decay  in 
the  neighbourhood  of  Q^.  For  simplicity  we  present  only  the  real  part  of  the  admissible  solutions 
of  (6.2)  although  it  is  clear  that  both  q2D  and  a  may,  in  general,  be  complex  while  q2D  is  always 
real.  Convergence  of  the  solution  q  towards  Qj  may  be  monitored  by  reference  to  either  the  local 
behaviour  of  the  solution  q  at  a  position  (sq,  yo)  on  S  or  by  monitoring  a  suitably  defined  global 
criterion  such  as  the  energy  contained  in  the  residuals  q2D;  alternatives  have  been  discussed  by 
Theofilis  (1998a).  Here  we  follow  the  first  approach  and  recover  the  parameter  a  by  monitoring 
the  solution  at  two  time-levels,  t  —  At  and  t,  where  At  may  but  need  not  be  the  time-step  in 
the  numerical  solution  algorithm.  Combining  (6.1)  and  (6.2)  it  follows  that  the  time-behaviour 
of  the  solution  may  be  monitored  by 

a  =  ln[qVq*-'^*]/At  «  dln[q*]/dt,  (6.3) 

where 
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The  approximation  in  (6.3)  holds  as  equality  in  the  case  of  linear  dependence  of  ln[q*]  on  time 
t.  Decay  of  residuals  is  indicated  by  cr  <  0.  A  first  statement  is  thus  in  place  without  reference 
to  a  particular  flow,  through  the  analytical  result  that  an  exponential  decay  of  residuals  near 
convergence  should  be  observed  as  a  consequence  of  the  separability  of  the  linearised  system  of 
equations  for  the  determination  of  residuals  in  time. 

(iii)  A  PHYSICAL  POINT  OF  VIEW  BASED  ON  NONPARALLEL  LINEAR  INSTABILITY  THEORY 

Explanation  of  the  further  observations  made  in  section  6  b  (i)  requires  calling  upon  an 
extension  of  the  classic  linear  instability  theory  proposed  by  Tollmien  (1929),  which  describes 
the  behaviour  of  small-amplitude  disturbances  superimposed  upon  an  one- dimensional  steady- 
state  basic  profile,  into  a  new  theory  which  is  concerned  with  small-amplitude  perturbations 
superimposed  upon  a  steady  two-dimensional  field.  In  so  doing,  the  many  and  often  questionable 
assumptions  related  with  the  so-called  parallel-flow  approximation  are  relaxed  and  the  linear 
instability  of  nonparallel  basic  states  may  be  analysed.  The  penalty  to  be  paid  in  resolving  two 
spatial  dimensions  simultaneously  is  the  need  for  numerical  solution  of  a  partial-derivative-based 
eigenvalue  problem  instead  of  the  straightforward  ordinary-differential-equation-based  system  of 
the  Orr-Sommerfeld  and  Squire  equations  (Drazin  and  Reid  1981).  One  of  the  early  successes 
of  the  nonparallel  two-dimensional  linear  instability  analysis  was  the  discovery  of  inviscid  short¬ 
wave  instability  of  two-dimensional  eddies  by  Pierrehumbert  (1986)  while  the  first  viscous  linear 
analysis  in  two  non-periodic  spatial  dimensions  known  to  us  is  the  work  of  Lee  et  al.  (1989)  on 
the  instability  of  flow  in  a  rectangular  enclosure  under  the  influence  of  gravity  and  temperature 
gradient.  More  recent  viscous  analyses,  in  step  with  modern  developments  in  algorithms  and 
hardware,  have  been  presented  by  Theofilis  (1998c). 

We  re-interpret  the  transient  solution  q  in  three-dimensional  physical  space  as  one  composed 
of  small-amplitude  three-dimensional  perturbations  q  =  {u,v,w,p)'^  superimposed  upon  Qj  = 
{u,v,w,p)'^,  the  latter  again  taken  to  be  two-dimensional.  Linearisation  about  Qj  is  permissible 
on  account  of  the  smallness  of  perturbations  compared  with  the  steady-state  and  the  resulting 
system  for  the  determination  of  q  is  separable  in  both  t  and  z  on  account  of  the  steadiness  and 
the  two-dimensionality  of  the  basic  flow  Q^.  Eigenmodes  are  introduced  in  these  directions  such 
that 

q{x,  y,  z,  t)  =  q{x,  y)  c.c.  (6.5) 

with  q  =  {u,v,w,p)'^  and  w  being  the  disturbance  velocity  component  in  the  z-direction. 
Complex  conjugation  is  introduced  in  (6.5)  since  q  is  real  while  all  three  of  q,  A  and  may  be 
complex.  In  the  framework  of  a  temporal  linear  nonparallel  instability  analysis  used  presently 
we  write  the  linearised  system  in  the  form  of  an  eigenvalue  problem  for  the  complex  quantity  fl, 
while  A  is  taken  to  be  a  real  wavenumber  parameter  describing  an  eigenmode  in  the  z-direction. 
The  real  part  of  is  related  with  the  frequency  of  the  instability  mode  while  its  imaginary  part 
is  the  growth/ damping  rate;  a  positive  value  of  fli  =  Im{fl}  indicates  exponential  growth  of  the 
instability  mode  q  in  time  t  while  fli  <  0  denotes  decay  of  q  in  time.  In  the  present  framework 
the  three-dimensional  space  comprises  S  extended  periodically  in  z  and  characterised  by  a 
wavelength  in  this  direction  which  is  associated  with  the  wavenumber  of  each  eigenmode.  A, 
through  Lz  =  27r/A. 

The  system  for  the  determination  of  and  q  takes  the  form  of  a  complex  nonsymmetric 
generalised  eigenvalue  problem 
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[£  -  {Vxu)]  u  -  {'Dyu)v  -  VxP  =  -ifi  u,  (6.6) 

-{Vxv)u  +  [C-  {Vyv)]  V  -  Vyp  =  -ifi  V,  (6.7) 

—  {'Dxw)u  —  {T>yw)v  +  Cw  —  iAp  =  — w,  (6-8) 

VxU  +  VyV  +  iXw  =  0  (6-9) 


subject  to  appropriate  boundary  conditions  on  dT,.  The  linear  operator 

£  =  (l/Re)  (Vl  +Vl-  A^)  -  uVx  -  vVy  -  \\w 

with  Vx  =  d/dx,  /dx"^,  Vy  =  d/dy  and  Vy  =  /dy^.  Comparison  of  (6. 1-6. 2)  and 

(6.5)  reveals  that  the  two  formalisms  are  related  in  the  limit  A  — 0.  However,  w  is  not  taken  a 
priori  to  vanish  within  the  framework  of  nonparallel  linear  instability;  three-dimensionality  of 
physical  space  is  preserved  and  the  existence  of  a  two-dimensional  steady-state  solution  Qj  is 
the  result  of  q  — 0  as  t  — oo.  The  comparison  of  (6. 1-6.2)  and  (6.5)  highlights  two  further  key 
ideas.  On  the  one  hand, 

residuals  acquire  the  physical  interpretation  of  one  of  the  linear  eigenmodes  which  pertain 

to  the  steady-state  Qj  and  have  A  =  0; 

on  the  other  hand, 

the  rate  of  decay  of  the  residuals  a  is  nothing  hut  the  damping  rate  fli  of  this  linear  pertur¬ 
bation, 

as  delivered  by  numerical  solution  of  the  partial- derivative  eigenvalue  problem  (6. 6-6. 9).  An¬ 
other  question  naturally  arising  concerns  the  physical  behaviour  of  the  system  when  the  least 
stable  member  of  the  linear  eigenspectrum  which  pertains  to  and  has  A  =  0  becomes  unstable. 
The  answer  is  clearly  that 

the  existence  of  an  unstable  (A  =  0)  — eigenmode  is  mutually  exclusive  with  the  ability  to 

obtain  a  converged  Q^. 

Prom  the  point  of  view  of  the  global  linear  instability  theory  based  on  the  partial  derivative 
eigenvalue  problem  (6. 6-6. 9)  the  unsteady  behaviour  of  two-dimensional  flow  may  be  related  to 
(A  =  0)— eigenmodes  approaching  conditions  of  neutral  stability  and  interacting  nonlinear ly. 

The  answer  to  the  second  question  posed  in  section  6  (a )  may  now  also  be  obtained  without 
reference  to  a  specific  flow  example.  The  existence  of  a  steady-state  Qj  in  a  2D  numerical 
simulation  is  synonymous  with  the  fact  that  all  (A  =  0)— eigenmodes  of  the  flow  have  Di  <  0. 
Modes  having  A  /  0,  on  the  other  hand,  may  be  either  growing  or  decaying  linearly.  In  case 
Di  <  0  V  A,  a  three-dimensional  numerical  simulation  performed  at  some  parameters  in  a  three- 
dimensional  domain  defined  by  S  and  an  arbitrary  periodic  extent  in  the  z-direction  will 
deliver  identical  results  for  a  converged  Qb  compared  with  that  of  a  two-dimensional  simulation 
performed  at  the  same  parameters  in  the  domain  S.  The  situation  changes  in  case  a  bracket  of 
wavenumbers  A  G  [Ai,A2]  exists  which  corresponds  to  unstable  modes.  The  largest  wavenumber 
A2  defines  a  length  £^2  =  27r/A2;  if  the  three-dimensional  simulation  is  performed  with  <  £^2 
again  no  difference  is  to  be  expected  between  its  result  for  Qj  and  that  of  a  two-dimensional 
simulation.  Both  will  converge  to  the  same  steady-state  solution  since  all  wavenumbers  of 
modes  defined  by  an  constrained  as  above  correspond  to  Di  <  0.  However,  if  >  £^2  at 
least  one  mode  in  the  three-dimensional  simulation  will  be  unstable,  which  will  result  in  the 
two-  and  three-dimensional  simulations  producing  different  solutions. 

We  return  to  the  observation  of  oscillatory  behaviour  of  the  residuals  near  convergence  and 
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differentiate  between  exponentially  decaying  residuals  of  either  sinusoidal  or  apparently  non¬ 
linear  nature.  A  linear  decay  of  ln[q*]  is  a  consequence  of  (A  =  0)— linear  eigenmodes  being 
stationary,  i.e.  having  fir  =  Re{fi}  =  0.  However,  other  stable  two-dimensional  member  of  the 
eigenspectrum  of  need  not  correspond  to  stationary  modes;  damped  travelling  modes  hav¬ 
ing  fir  /  0  will  manifest  themselves  in  the  time-accurate  simulation  as  residuals  of  sinusoidal 
character  the  magnitude  of  which  decays  exponentially.  On  the  other  hand,  the  unambiguously 
linear  dependence  of  ln[q*]  on  t  in  the  neighbourhood  of  Qj  is  the  consequence  of  the  existence  of 
a  spectrum  comprising  modes  which  are  clearly  separated  in  parameter  space  from  one  another. 
The  co-existence  of  several  two-dimensional  eigenmodes  of  approximately  the  same  damping  rate 
can  lead  to  their  nonlinear  interaction  and  difficulty  to  observe  a  behaviour  governed  by  non¬ 
parallel  linear  instability  theory.  Comparison  of  power  spectral  analysis  of  the  time-dependent 
DNS  signal  and  the  results  of  the  partial- derivative  eigenvalue  problem  (6. 6-6. 9)  may  shed  light 
upon  the  two-dimensional  eigenmodes  involved  in  such  a  nonlinear  interaction. 

(iv)  On  the  time  oe  integration  until  convergence 

Straightforward  rearrangement  of  (6. 1-6.2)  delivers  an  estimate  of  the  time  necessary  (under 
linear  conditions)  for  the  least  stable  global  mode  present  in  the  numerical  solution  to  be  reduced 
from  an  amplitude  Aq  to  a  lower  level  Ai,  which  may  be  calculated  from 

Ta,/a,  =  ln(Ai/Ao)/(-fii),  (6.10) 

where  fii  is  the  damping  rate  of  the  mode  in  question.  The  worst  case  scenario  in  a  time- 
accurate  integration  is  that  the  solution  will  be  attracted  by  the  least-stable  global  eigenmode 
developing  upon  Qj  and  having  A  =  0  throughout  the  course  of  the  simulation.  An  upper  bound 
for  the  time  necessary  for  the  steady-state  to  be  obtained  may  then  be  offered  by  (6.10)  in  which 
fii  is  the  damping  rate  of  this  mode.  Defining,  for  example,  convergence  as  the  reduction  of  an 
0(1)  residual  by  10  orders  of  magnitude  results  in  an  integration  time  of  T^q-io  «  23/|fii|.  This  is 
a  conservative  estimate  since  it  is  occasionally  observed  that  other  stronger  damped  eigenmodes 
will  come  into  play  early  in  the  simulation  and  the  least-damped  eigenmode  will  only  determine 
the  late  stages  of  the  convergence  process. 

An  associated  point  concerns  the  misconception  which  often  exists  that  initialising  the  numer¬ 
ical  solution  for  Qj  at  some  Reynolds  number  from  a  state  which  is  ’close’  to  the  one  desired, 
for  instance  using  the  converged  solution  at  a  somewhat  different  Reynolds  number,  may  reduce 
the  integration  time  in  the  context  of  a  time-accurate  solution.  The  present  analysis  shows  this 
to  be  a  misplaced  expectation.  If  there  exists  an  0(1)  deviation  between  the  target  solution  and 
its  initial  estimate,  the  deviation  has  to  be  reduced  in  magnitude  during  an  integration  of  the 
equations  of  motion  for  the  length  of  time  determined  by  the  least  damped  two-dimensional 
(A  =  0)— eigenmode  at  the  specific  Reynolds  number.  It  is  this  eigenmode  of  the  flow  and  not 
the  initial  state  which  determines  the  length  of  the  integration  time  for  Qj.  The  ideas  dis¬ 
cussed  in  the  previous  subsection  on  the  other  hand,  lead  to  an  algorithm  application  of  which 
may  save  substantial  amounts  of  the  integration  time  necessary  for  reduction  of  residuals  to 
machine-roundoff  level. 

(v)  Recovery  oe  the  converged  solution  erom  transient  data 

Having  identified  small-amplitude  residuals  in  the  calculation  as  the  least  damped  global  two- 
dimensional  eigenmodes  of  the  flow,  it  is  now  possible  to  utilise  this  information  in  order  to 
recover  the  converged  steady-state  solution  from  transient  data,  without  having  to  pursue  the 
time  integration  of  the  equations  of  motion  until  convergence  in  time  is  obtained.  Combining 
(6.1),  (6.2)  and  (6.5)  one  obtains 

q{x,y,t)  =  Qb{x,y)  +  e  j^qr  cos  fir<  —  qi  sin  firt j  e*’^*,  (6-11) 
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where  %  =  Re{q},  qi  =  Im{q}  and  q  is  one  of  the  (A  =  0)— linear  eigenmodes  in  (6.5).  It 
should  be  stressed  here  that  the  following  discussion  is  applicable  to  transient  data  for  which 
(6.11)  holds,  namely,  solutions  for  which  the  entire  time-dependence  of  the  solution  is  exhibited 
in  the  residuals;  in  other  words,  the  present  analysis  is  based  on  the  self-consistent  premises 
that  dQijfdt  =  0.  Further,  it  is  noted  that  q  may  but  need  not  be  the  least-damped  member 
of  the  eigenspectrum  of  Qj;  the  only  prerequisite  for  the  validity  of  the  following  discussion  is 
that  the  transient  solution  has  reached  a  regime  of  exponential  decay  of  residuals.  A  final  point 
is  that  the  signal  near  convergence  need  not  be  composed  of  a  single  damped  eigenmode  as 
(6.11)  implies.  However,  the  elements  of  the  theory  for  the  recovery  of  from  a  signal  being 
composed  of  several  stationary  (fir  =  0)  and  travelling  (fir  /  0)  linearly  damped  eigenmodes 
may  be  exposed  by  reference  to  (6.11)  on  which  we  focus  our  attention. 

The  calculation  of  Qj  from  transient  data  for  q  follows  in  two  stages.  First,  elementary  signal 
analysis  techniques  deliver  the  results  for  fir  and  a.  Second,  once  fir  and  a  have  converged  in 
time  (6.11)  may  be  used  to  calculate  Q^.  The  circular  frequency  fir  is  calculated  from  the  the 
period  of  oscillations  in  the  time-signal  of  q  which,  in  turn,  is  identified  by  the  maxima  in  the 
signal.  Independently,  in  order  to  calculate  a  we  re-write  (6.11)  as 


W 


dq  ^  d^q 


+  (a^+fir^)^-2a 


=  0. 


(6.12) 


This  expression  may  be  evaluated  at  those  times  that  dq/dt  =  0  in  the  course  of  the  time- 
integration,  i.e.  at  the  same  times  that  fir  is  calculated.  At  these  times  the  magnitude  of  a  is 
given  by 


a 


1  {d^q/dt^) 


2  (d^q/dt^) 


{dq/dt)=0 


(6.13) 


In  case  fir  =  0,  a  monotonic  dependence  of  dqfdt  on  t  is  usually  observed  from  the  beginning 
of  the  calculation  until  convergence,  with  dqfdt  =  0  only  at  convergence.  In  this  case,  the 
magnitude  of  a  may  be  calculated  using 


{d\ldt^) 

{dqfdt) 


(6.14) 


With  (T  and  fir  converged  in  time  (6.11)  may  be  written  as  a  linear  system  of  three  equations 
at  three  times  ti  =  t,t2  =  t  +  At  and  tz  =  t  +  2 At  for  three  unknowns,  Qj,  qr  and  qi  with 
the  transient  solution  q„  =  q{x,y,tn)  known  at  these  times.  Simple  algebra  delivers  the  desired 
converged  steady-state  solution  as 


Qb 


qi  -  2  q2  e*^^*  cos  fir  At  qz 
g2<TAt  _  2  cos  fir  At  -|-  1 


(6.15) 


As  an  aside,  the  spatial  structure  (qrjqi)  of  the  linear  eigenmode  q  may  also  be  recovered  to 
within  an  arbitrary  constant  from  the  same  linear  system.  Equivalently,  if  only  the  converged 
steady-state  solution  is  of  interest,  the  expression 


may  be  used  for  the  recovery  of  Qj  from  transient  data  for  q  and  its  first  two  time-derivatives. 
Either  of  (6.15)  or  (6.16)  may  be  used  for  the  cases  of  residuals  corresponding  to  stationary 
(fir  =  0)  or  travelling  (fir  /  0)  single  linear  eigenmodes. 

This  idea  may  be  extended  to  extract  Qj  from  a  DNS  signal  comprising  several  linearly 
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decaying  eigenmodes  superimposed  upon  the  steady-state  solution, 

q  =  Q6  +  Y^£n{qn,rCOsQ.n,rt  “  sin  (6.17) 

n 

As  an  example,  in  the  case  of  one  stationary 

eiqi.re””'*  (6.18) 

and  one  travelling 

£2(q2,r  cos  firt  —  q2,i  sin  (6.19) 

linear  disturbance  being  present  in  the  signal,  one  may  first  extract  information  for  the  damp¬ 
ing  rate  of  the  stationary  mode  from  the  signal  itself  and  for  the  damping  rate  and  frequency  of 
the  travelling  disturbance  from  the  first  time-derivative  of  the  DNS  signal  for  q.  Subsequently, 
one  may  solve  the  {2NxNy)  x  {2NxNy)  system  defined  by  writing 

((t|  -|-  flr)Q6  +  (^1  +  cr|  -|-  fly  “  2(TiO'2)eiqi,re‘^^*  = 

“  2(72-^ -I- ((t| -I- flr)q  (6.20) 

at  two  consecutive  times  ti  and  t2  for  Q^,  and  eiqi,^?  where  Nx  and  Ny  are  the  number  of 
points  discretising  the  x—  and  y— spatial  directions,  respectively. 

The  accuracy  by  which  fir  and  a  are  determined  depends  on  that  by  which  the  first  three 
time-derivatives  of  q  are  calculated;  this,  in  turn,  depends  on  the  time-step  in  the  calculation 
and  the  number  of  fields  stored  in  order  for  backward  differentiation  formulae  to  be  applied. 
Since  the  time-step  is  controlled  by  CFL  considerations,  it  is  advisable  to  store  a  reasonably 
high  number  of  fields  in  order  for  high  accuracy  of  fir  and  a  and,  in  turn,  of  Qj  to  be  obtained. 
The  calculations  to  be  presented  in  what  follows  have  been  performed  using  five-point  backward 
differencing  formulae  on  an  equidistant  grid  (Abramowitz  and  Stegun  1956). 

At  conditions  at  which  a  steady-state  solution  exists  most  two-dimensional  global  eigenmodes 
of  the  converged  steady-state  are  heavily  damped  ((t„  =  0(1)  in  equation  (6.17)).  Consequently, 
if  the  time-integration  of  the  equations  of  motion  is  pursued  long  enough,  only  a  handful  of 
(A  =  0)— global  eigenmodes  will  survive  and  persist  in  the  DNS  signal.  Clearly,  it  is  the  least 
damped  of  the  global  instabilities  that  will  determine  the  ultimate  behaviour  of  the  solution.  In 
determining  whether  one  integrates  the  equations  of  motion  until  all  but  the  least-damped  of 
the  eigenmodes  have  subsided  in  order  to  apply  (6.15)  or  (6.16)  or  one  recovers  Qj  at  an  earlier 
time  from  a  signal  in  which  a  number  of  damped  eigenmodes  still  persist  one  should  take  into 
account  the  following  factors. 

First,  the  efficiency  of  the  specific  DNS  algorithm  determines  whether  the  cost  of  integrating 
the  equations  of  motion  until  convergence  is  acceptable  at  given  flow  parameters.  The  cost  of 
computing  fir,  cr,  intermediate  values  of  and  monitoring  convergence  of  all  these  quantities, 
possibly  for  several  eigenmodes,  must  also  be  weighed  against  the  straightforward  approach  of 
pursuing  the  time-integration  in  the  DNS  until  convergence.  However,  at  all  Reynolds  numbers 
studied  in  the  prototype  flow  monitored  both  a  and  fir  of  individual  modes  have  converged  within 
the  first  quarter  to  half  of  the  total  integration  time,  making  further  time- integration  superfluous. 
While  the  integration  time  until  convergence  is  short  at  low  Reynolds  numbers,  on  account 
of  large  damping  rates  of  the  least-damped  linear  eigenmodes,  at  increasingly  large  Reynolds 
numbers  the  magnitude  of  the  damping  rates  becomes  increasingly  smaller  and  application  of 
the  ideas  exposed  in  this  section  becomes  increasingly  attractive  in  order  for  substantial  savings 
in  computing  effort  to  be  materialised. 
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(c)  Results  for  the  square  lid-driven  cavity 

An  example  flow  in  which  these  ideas  may  be  illustrated  is  the  classic  lid-driven  cavity 
(Burggraf  1966).  In  its  function  as  a  testbed  for  numerous  algorithms  this  flow  has  generated  a 
substantial  amount  of  information  which  is  relevant  to  the  preceding  discussion.  Calculations  for 
Qb  were  performed  using  a  two-dimensional  spectrally-accurate  algorithm  for  direct  numerical 
simulation  of  flow  in  nonperiodic  geometries.  The  code  is  based  on  a  real-space  eigenvalue- 
decomposition  of  the  spectral  collocation  differentiation  matrices  extending  ideas  discussed  by 
Ku  et  al.  (1988)  and  uses  one  member  of  the  low-storage  second-order  accurate  time-integration 
schemes  put  forward  by  Spalart  et  al.  (1991).  A  spectral  algorithm  was  chosen  in  order  for  opti¬ 
mal  accuracy  to  be  obtained  on  a  low  number  of  collocation  points,  the  latter  being  dictated  by 
the  maximum  number  of  points  on  which  numerical  solution  of  the  partial-derivative  eigenvalue 
problem  is  feasible  using  current  computer  technology.  Solutions  were  obtained  using  Jacobi 
polynomials  for  the  spatial  discretisation  at  resolutions  depending  on  the  Reynolds  number  and 
ranging  from  32^  to  128^  spectral  collocation  points.  The  time-steps  at  the  different  Reynolds 
numbers  were  kept  well  below  those  dictated  by  the  CFL  condition  in  order  for  reasonable  ac¬ 
curacy  of  the  results  of  cr  to  be  ensured.  In  view  of  our  arguments  being  based  on  nonparallel 
linear  instability  analysis  and  the  well-known  sensitivity  of  linear  instability  analysis  results  on 
the  accuracy  of  the  basic  flow,  we  first  present  a  validation  of  both  the  basic  flow  and  the  partial 
derivative  eigenvalue  problem. 

(i)  Validation  studies 

The  accuracy  of  the  converged  steady-state  solutions  is  first  assessed  by  comparison  with 
the  established  works  of  Ghia  et  al.  (1982)  and  Schreiber  and  Keller  (1983).  Converged  basic 
states  have  been  calculated  at  several  Reynolds  numbers  of  which  we  present  calculations  at 
Re  =  400, 1000, 3200  and  4000,  the  first  three  obtained  on  32^  and  the  last  on  48^  Legendre 
collocation  points.  At  Re  =  400  and  1000  both  aforementioned  works  present  results  while 
at  the  higher  Reynolds  numbers  we  compare  our  calculations  individually  with  either  work. 
Interestingly,  aside  from  the  locations  and  maximum  values  of  streamfunction  and  vorticity  in 
the  primary  vortex  core,  Schreiber  and  Keller  (1983)  analysed  and  presented  their  results  in  the 
form  of  a  converging  series  calculated  by  Richardson  extrapolation.  Comparisons  are  presented 
in  a  twofold  manner.  The  comparison  of  our  calculations  for  the  location  and  maxima  in  the 
stream-function  ip  and  the  vorticity  ^  with  those  of  the  reference  works  is  shown  in  Table  3; 
note  that  Schreiber  and  Keller  (1983)  define  ^  to  have  an  opposite  sign  to  that  of  Ghia  et  al. 
(1982)  and  the  present  work.  Although  the  overall  agreement  of  all  results  is  quite  reasonable 
marginal  differences  exist.  These  may  be  attributed  to  the  different  grids  used  in  all  three  works, 
making  an  interpolation  procedure  necessary  for  detailed  comparisons.  To  this  end,  we  employed 
a  piecewise  cubic  procedure  to  transfer  our  results  onto  the  (different)  maxima  of  the  benchmark 
calculations.  Our  interpolated  values  as  well  as  the  results  of  Ghia  et  al.  (1982)  and  Schreiber 
and  Keller  (1983)  are  presented  in  Table  4b,  where  the  individual  comparisons  demonstrate  a 
substantially  more  satisfactory  agreement  of  our  calculations  with  both  benchmark  works  at 
low  i?e— values  and  especially  with  the  Richardson  extrapolated  results  of  Schreiber  and  Keller 
(1983)  at  the  highest  Reynolds  number  monitored. 

It  is  well-known  from  comparisons  of  three-dimensional  DNS  results  and  one-dimensional 
Orr-Sommerfeld-based  linear  instability  analysis  that  details  of  the  steady  basic  state  strongly 
influence  the  accuracy  of  the  growth/ damping  rates  of  linear  eigenmodes  (Kleiser  and  Schumann 
1980).  The  remaining  differences  between  our  results  for  the  two-dimensional  steady-states  in 
the  lid-driven  cavity  and  those  of  the  benchmark  works  are  next  assessed  in  this  light,  from  the 
point  of  view  of  their  influence  on  the  global  linear  instability  analysis  results.  Two  solutions 
of  the  partial  derivative  eigenvalue  problem  (6. 6-6. 9)  for  the  lid-driven  cavity  exist,  those  of 
Ramanan  and  Homsy  (1994)  (RH)  and  Ding  and  Kawahara  (1998)  (DK).  These  authors  have 
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presented  linear  instability  analyses  of  the  square  lid-driven  cavity  flow  which  deliver  consistent 
results  at  low  Re  but  predict  different  critical  Reynolds  number  values  for  linear  amplification 
of  three-dimensional  perturbations.  While  individual  comparisons  are  certainly  possible,  at  high 
Reynolds  numbers  neither  work  presents  results  for  the  two-dimensional  global  linear  instabilities 
which  are  central  to  the  theme  of  the  present  paper.  We  therefore  refrain  here  from  discussion 
of  three-dimensional  linear  instability  and  the  issue  of  a  linear  critical  Reynolds  number  and 
monitor  a  low  value  of  the  Reynolds  number,  Re  =  200,  at  which  both  RH  and  DK  present 
results  at  A  =  0. 

Table  5  shows  the  tabulated  values  of  RH,  the  graphically  reproduced  results  of  DK  and  our 
solutions  of  the  partial-derivative  eigenvalue  problem  (6. 6-6. 9).  The  overall  agreement  of  the 
previous  and  the  present  instability  analyses  is  quite  good  and  all  results  indicate  the  experi¬ 
mentally  established  fact  of  stability  of  the  two-dimensional  flow  in  the  lid-driven  cavity  at  this 
Reynolds  number  (Burggraf  1966).  Regarding  the  quality  of  the  basic  flow,  it  may  be  inferred 
from  the  results  of  Table  5  that  the  basic  states  of  both  RH  and  DK  and  the  present  work  are 
practically  identical  for  the  purposes  of  the  linear  instability  analysis  that  follows,  at  least  at 
the  Reynolds  numbers  discussed. 


(ii)  Numerical  residuals  and  (A  =  0)  linear  eigenmodes  in  the  square  lid-driven 

CAVITY 

Figs.  26-30  show  the  convergence  histories  of  the  two-dimensional  DNS  at  several  Reynolds 
numbers,  with  the  qualitative  behaviour  of  residuals  discussed  earlier  observed.  The  convergence 
of  the  rate  of  decay  of  residuals  a,  calculated  using  (6. 3-6.4),  is  shown  in  Table  6  at  Re  =  100, 200 
and  300.  Also  shown  is  the  damping  rate  fli  of  the  least-damped  eigenmode  having  A  =  0  as 
obtained  by  linear  analysis,  based  on  the  partial-derivative  eigenvalue  problem  (6. 6-6. 9),  of  the 
converged  steady-state  Qb  corresponding  to  each  Reynolds  number.  The  excellent  agreement 
between  the  two  quantities  leaves  little  room  for  doubt  that  numerical  residuals  may  be  identified 
as  being  the  least-damped  (A  =  0)— eigenmode  of  the  corresponding  converged  steady-state.  It 
is  interesting  to  note  here  that  such  an  agreement  could  not  be  obtained  when  we  followed  the 
commonly-used  procedure  to  terminate  the  steady-state  calculation  after  a  decay  of  residuals  by 
an  arbitrarily  defined  seemingly  adequate  small  number  of  orders  of  magnitude,  say  5-6.  Such 
poorly  converged  in  time  basic  states  may  be  viewed  as  comprising  a  small  unsteady  component 
the  linear  instability  analysis  of  which  is  bound  to  deliver  erroneous  results.  Further,  it  is  worth 
mentioning  that  the  prediction  (6.10)  of  the  time  necessary  to  integrate  the  equations  of  motion 
until  convergence  in  time  is  in  line  with  the  results  of  Table  6  and  Fig.  26. 

A  clearly  defined  single  value  of  a  which  determines  the  behaviour  of  residuals  in  the  entire 
course  of  the  time-integration  is  a  result  of  a  two-dimensional  eigenspectrum  of  Qb  in  which 
the  least  damped  two-dimensional  (A  =  0)— eigenmodes  are  stationary  and  well  separated  in 
parameter  space  from  their  more  stable  counterparts.  The  situation  becomes  more  intricate, 
but  still  amenable  to  analysis,  as  the  Reynolds  number  increases.  Qualitative  differences  may 
be  found  between  the  results  of  Figs.  26  and  27,  although  all  simulations  were  started  from 
the  same  initial  condition  ip  =  (^  =  0.  While  in  both  sets  of  results  a  short  initial  transient  is 
followed  by  exponential  decay  of  residuals,  in  the  first  set  this  decay  pursues  at  the  same  rate 
for  almost  two  decades  while  in  the  second  two  different  rates  of  exponential  decay  of  residuals 
are  demonstrated.  Inspection  of  the  full  spectra  delivered  by  numerical  solution  of  (6. 6-6. 9) 
at  each  Reynolds  number  reveals  that  as  the  Reynolds  number  increases  an  increasingly  larger 
number  of  eigenmodes,  both  stationary  and  travelling  appear  in  the  eigenspectrum  of  Qb,  having 
damping  rates  approximately  equal  with  that  of  the  least  damped  eigenmode.  As  a  consequence, 
the  numerical  solution  may  initially  be  attracted  to  a  different  than  the  least  damped  (A  = 
0)— eigenmode  but  its  long-time  behaviour  will  be  determined  by  the  latter  disturbance.  In  both 
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the  Re  =  500  and  the  Re  =  1000  results  of  Fig.  27  the  damping  rates  of  the  least  and  the  next 
more  stable  mode  are  presented  as  symbols  superimposed  upon  the  curves  used  to  determine  a. 

Yet  another  qualitatively  different  behaviour  is  observed  in  the  time-signal  of  q  as  a  con¬ 
sequence  of  a  further  increase  of  the  Reynolds  number.  Alongside  the  least  damped  stationary 
mode  travelling  disturbances  appear,  as  seen  in  the  results  of  Figs.  28-30.  In  all  three  figures  a{t) 
assumes  the  form  of  exponentially  decaying  disturbances.  However,  while  at  the  lowest  Reynolds 
number  a  clearly  identifiable  sinusoidal  perturbation  may  be  seen,  having  fir  ~  0.97  ±  0.01,  a 
barely  perceptible  deviation  from  a  single  oscillatory  disturbance  (fir  ~  0.954  ±  0.012)  may  be 
seen  at  Re  =  5000;  at  Re  =  7500  the  solution  demonstrates  a  behaviour  which  might  be  inter¬ 
preted  either  as  nonlinearity  or  as  superposition  of  two  exponentially  decaying  linear  sinusoidal 
disturbances  having  frequencies  of  fir  =  0.933  and  0.945.  In  order  to  analyse  these  observations 
we  pursue  two  independent  paths.  First,  we  perform  a  nonparallel  linear  instability  analysis  of 
the  converged  steady  state  at  each  Reynolds  number  and  monitor  the  least-stable  member  of 
the  eigenspectrum,  which  turns  out  to  be  a  stationary  linear  eigenmode.  Second,  we  perform 
a  discrete  Fourier  transform  (DFT)  of  the  DNS  signal  for  q  at  Re  =  2500, 5000  and  7500  and 
compare  the  results  with  the  eigenvalues  of  the  travelling  disturbances  delivered  by  the  linear 
instability  analysis. 

Table  7  shows  that  a  progressive  deviation  of  the  rate  of  decay  of  the  residuals  from  the 
damping  rate  of  the  least-damped  eigenmode  occurs  as  Re  increases.  This  result  suggests  that 
as  the  Reynolds  number  increases  nonlinear  interaction  of  the  least  stable  eigenmodes  may  cause 
a  departure  of  the  numerical  solution  from  a  behaviour  predicted  by  nonparallel  linear  theory. 
The  role  that  the  least  stable  members  of  the  full  eigenvalue  spectrum  play  in  the  dynamics  of 
the  flow  may  be  inferred  from  the  results  of  Figs.  31-  33.  In  Fig.  31  we  present  the  DFT  of  the 
DNS  signal  for  ■0(0.5, 0.5)  scaled  by  the  maximum  value  of  the  spectral  density.  A  single  peak  at 
27r/  «  1,  albeit  of  somewhat  wide  support,  dominates  over  two  much  smaller  peaks  at  27r/  =  0 
and  27r/  «  2.  Shown  are  also  the  results  of  (6. 6-6. 9)  for  fir,  arbitrarily  placed  on  the  vertical 
axis  for  readability.  An  one-to-one  correspondence  between  the  peaks  in  the  spectrum  and  the 
values  of  fir  for  stationary  and  travelling  linear  eigenmodes  is  clearly  identifiable.  Interestingly, 
the  width  of  the  support  of  the  peaks  is  found  to  be  associated  with  the  existence  of  more 
than  one  eigenvalues  in  the  partial-derivative  eigenvalue  problem  spectrum,  at  both  fir  «  1  and 
fir  ~  2.  The  origin  of  the  existence  of  only  harmonics  of  the  first  travelling  eigenmode  in  the 
full  eigenvalue  spectrum  deserves  further  investigation.  It  should  be  noted  here  that  a  Krylov 
subspace  iteration  method  has  been  used  for  the  solution  of  the  partial- derivative  eigenvalue 
problem,  which  results  in  only  a  window  of  the  eigenvalue  spectrum  being  captured  at  any 
single  calculation.  The  number  of  converged  eigenvalues  recovered  increases  as  the  subspace 
dimension  increases.  However,  the  neighbourhood  of  (fir,  fii)  =  0  has  been  well  resolved  in  all 
results  presented  here.  A  higher  Krylov  subspace  dimension  has  been  found  to  deliver  additional 
eigenvalues  at  higher  frequencies. 

While  the  agreement  between  the  frequencies  in  the  DNS  signal  and  those  of  the  nonparallel 
linear  analysis  of  the  converged  steady  state  is  evident,  the  results  of  Fig.  31  do  not  provide 
any  information  on  the  damping  rates  fii  of  the  disturbances  whose  frequency  lies  at  fir  ~  1  in 
relation  to  those  at  different  frequencies.  For  our  argument  that  the  residuals  in  the  calculation 
may  be  identified  as  the  least  stable  of  the  two-dimensional  (A  =  0)— eigenmodes  to  be  valid,  the 
damping  rate  of  the  linear  disturbances  at  fir  ~  1  must  be  lower  than  that  of  modes  with  higher 
frequencies;  this  is  a  point  to  which  we  will  return  shortly.  Qualitatively  analogous  results  are 
obtained  at  Re  =  5000,  seen  in  Fig.  32.  Besides  the  slight  shift  towards  lower  frequencies,  the 
quantitative  difference  with  the  results  at  Re  =  2500  is  that  the  strength  of  the  eigenmodes  at 
fir  ~  2  is  substantially  larger  than  that  of  their  counterparts  at  Re  =  2500  in  relation  to  the 
strength  of  the  respective  modes  at  fir  ~  1;  this  is  the  origin  of  the  slight  deviation  from  a  purely 
sinusoidal  behaviour  of  the  signal  at  this  Reynolds  number,  seen  in  Fig.  29.  Further,  additional 
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eigenmodes  appear  alongside  the  counterparts  of  those  seen  at  Re  =  2500  at  fir  =  0  and  fir  «  1 
and  new  modes  appear  at  fir  ~  3  and  fir  ~  4.  Finally,  at  Re  =  7500,  the  pattern  discovered  at 
the  lower  Reynolds  number  values  qualitatively  repeats  itself,  with  the  appearance  of  additional 
modes  at  the  same  and  new  at  high-frequencies;  furthermore  a  new  mode  which  does  not  fit  in 
the  period-doubling  scenario  discussed  also  is  present  in  the  linear  instability  analysis  results, 
which  the  DFT  reveals  to  be  too  weak  to  play  an  important  role  in  the  dynamics  of  the  flow  at 
this  Reynolds  number  value. 

We  return  to  the  question  of  damping  rates  of  the  linear  instability  modes  and  present  in 
Fig.  34  the  full  spectrum  of  eigenvalues  in  the  neighbourhood  of  fi  =  0  at  Re  =  2500, 5000,  and  7500. 
Results  of  signiflcance  in  this  figure  are  the  following.  First,  as  the  Reynolds  number  increases 
the  flow  becomes  less  stable  to  two-dimensional  linear  (A  =  0)— eigenmodes.  Second,  in  all  three 
Reynolds  numbers  the  least  stable  modes  are  stationary  disturbances.  Third,  perfect  symmetry 
about  fir  =  0  may  be  observed  in  the  results,  as  should  be  expected  from  the  ability  to  refor¬ 
mulate  (6. 6-6. 9)  as  a  real  eigenvalue  problem.  Consistent  with  the  DFT  results  of  the  signal 
discussed  earlier,  the  eigenmodes  at  fir  ~  1  are  less  stable  than  their  counterparts  at  higher  fre¬ 
quencies.  Comparing,  for  example,  the  Re  =  5000  eigenmodes  (fir,fii)  =  (0.967,-0.0158)  and 
(fir,  fii)  =  (1.921,  —0.0319)  one  finds  that,  if  introduced  at  the  same  initial  amplitude  in  the  flow, 
the  second  mode  would  be  reduced  by  a  given  number  of  orders  of  magnitude  in  amplitude  in 
approximately  half  as  long  an  integration  time  as  that  required  for  the  first  mode  to  experience 
the  same  reduction  of  amplitude. 

(iii)  The  critical  Reynolds  number  oe  (A  =  0)- linear  disturbances 

The  preceding  discussion  leads  to  re-examination  of  the  question  of  a  critical  Reynolds  number 
for  linear  growth  of  two-dimensional  global  instabilities  in  the  square  lid-driven  cavity.  Consistent 
with  well-established  numerical  solutions  for  the  steady-state  in  this  flow  the  nonparallel  linear 
instability  analysis  results  of  §  6  (c )  ii  deliver  a  least  damped  stationary  (A  =  0)— eigenmode 
which  has  a  damping  rate  whose  magnitude  decreases  with  increasing  Reynolds  number.  The 
dependence  of  fii  on  Re  for  this  mode  has  been  obtained  at  several  Reynolds  numbers  and  is 
presented  by  symbols  in  Fig.  35.  Analysis  of  the  results  for  the  damping  rate  of  the  least-damped 
eigenmode  as  function  of  the  Reynolds  number  delivers  a  curve-fit  of  the  data  by  using 

fii  =  -109.071  (6.21) 

The  curve  defined  by  (6.21)  is  also  shown  in  Fig.  35  by  a  solid  line  and  has  been  found  to 
deliver  reasonably  accurate  predictions  of  fii  at  Re  >  1000,  where  the  calculated  data  may  be 
collapsed  onto  a  single  curve.  The  upper  bound  of  the  Reynolds-number  range  in  which  (6.21) 
may  be  used  with  confidence  to  predict  the  rate  of  decay  of  residuals  and  the  associated  time  of 
integration  of  the  equations  of  motion  until  a  steady-state  solution  is  reached  must  be  Re  «  10^, 
a  value  below  which  a  multitude  of  two-dimensional  numerical  solutions  have  demonstrated  the 
existence  of  converged  two-dimensional  states.  In  the  framework  of  the  current  nonparallel  linear 
instability  analysis  this  should  manifest  itself  by  Qj  losing  its  stability  in  a  linear  framework 
to  amplified  two-dimensional  perturbations  having  A  =  0.  However,  as  has  been  mentioned,  the 
existence  of  a  converged  steady-state  solution  is  synonymous  with  all  global  eigenmodes  of  the 
flow  being  stable.  Another  possibility  is  that  the  nonlinear  interaction  of  two-dimensional  global 
neutrally-stable  disturbances  as  the  Reynolds  number  increases  may  be  held  responsible  for  the 
observed  inability  to  obtain  a  converged  steady  state  solution.  However,  from  (6.21)  it  follows 
that  fii  <  0,  VRe  and  the  flow  remains  stable  to  all  two-dimensional  (A  =  0)— eigenmodes.  Two 
aspects  of  this  prediction  should  be  stressed  here.  First,  (6.21)  is  a  curve-fit,  at  best  valid  up 
to  the  highest  Reynolds  number  used  to  produce  it.  Re  =  7500.  Second,  the  filling-up  of  the 
eigenspectrum  and  the  associated  nonlinear  interaction  of  some  of  the  least  stable  eigenmodes 
as  the  Reynolds  number  increases,  causes  a  systematic  departure  of  the  numerical  solution  for  q 
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from  one  determined  by  a  single  eigenmode  of  the  nonparallel  linear  instability  theory,  as  already 
shown  in  the  results  of  Table  7.  On  the  other  hand,  the  trend  predicted  by  (6.21)  is  correct, 
namely  that  the  damping  rates  of  two-dimensional  global  linear  instabilities  as  Re  increases 
are  exponentially  small  in  magnitude.  As  such,  an  increasingly  large  number  of  global  modes 
may  be  considered  neutrally  stable  at  large  Reynolds  numbers;  it  is  therefore  likely  that  the 
second  scenario,  namely  nonlinear  interaction  of  near  neutrally  stable  two-dimensional  global 
flow  eigenmodes  is  responsible  for  the  observed  loss  of  ability  to  obtain  a  steady-state  solution 
of  the  equations  of  motion  at  Re  >  10‘*  (f.e.  E  and  Liu  1996). 

(iv)  Obtaining  the  converged  steady- state  solution  erom  non- converged 

TRANSIENT  DATA 

The  preceding  discussion  has  demonstrated  the  association  of  residuals  in  two-dimensional 
incompressible  DNS  calculations  with  the  two-dimensional  global  linear  instability  modes  of  the 
converged  steady  state.  In  this  section  we  present  examples  of  recovery  of  steady-state  solutions 
from  transient  DNS  data  using  this  information  and  the  algorithm  of  §  6  (b )  v.  We  stress  that 
the  applicability  of  the  algorithm  is  intimately  linked  with  the  quality  of  the  DNS  and  the  initial 
conditions  used  for  the  simulation,  since  both  determine  when,  for  what  length  of  time  and  to 
which  linear  eigenmode  the  time-accurate  solution  will  be  attracted  in  the  course  of  the  time- 
integration.  Here  we  present  a  discussion  of  some  parameters  which  affect  the  results  returned 
by  the  algorithm  in  a  few  Reynolds  number  cases  of  those  on  which  the  algorithm  was  validated. 

Results  at  Re  =  100  and  1000  are  shown  in  Table  8;  at  each  Reynolds  number  we  have 
performed  three  sets  of  calculations,  two  direct  numerical  simulations  and  one  solution  of  the 
partial-derivative  eigenvalue  problem.  Both  DNS  start  from  the  initial  condition  ip  =  (^  =  0 
for  the  flow  streamfunction  and  vorticity,  respectively.  On  the  one  hand,  the  converged  ’exact’ 
steady  state  Qj  has  been  calculated  by  marching  the  equations  of  motion  until  such  a  time  i 
that  the  residuals  were  reduced  to  machine- roundoff  level,  using  64-bit  arithmetic  and  monitoring 
convergence  along  the  lines  discussed  in  §  6  (c )  i.  On  the  other  hand,  we  have  run  another  DNS 
but  marched  the  equations  of  motion  until  such  a  time  t  was  reached  at  which  a  linear  regime 
was  identifled  by  the  convergence  in  time  of  fir  (when  applicable)  and  a.  The  time-marching 
was  then  interrupted  and  either  (6.15)  or  (6.16)  was  solved  for  the  respective  ’estimated’  steady- 
state  solution  q.  Finally,  the  partial-derivative  eigenvalue  problem  (6. 6-6. 9)  was  solved  for  two- 
dimensional  disturbances  (A  =  0)  developing  upon  Qj  and  the  eigenvalue  spectrum  pertaining  to 
the  flow  at  each  Reynolds  number  was  recovered.  The  results  were  compared  both  in  terms  of  the 
magnitude  of  the  relative  discrepancy  of  the  two  DNS-obtained  solutions  Aq  =  |(q  —  Q6)/Q6| 
and  by  monitoring  the  difference  between  a  in  the  second  set  of  DNS  and  fii.  Table  8  shows  the 
resolutions  and  time-steps  used  in  several  simulations,  the  time  t  at  which  a  converged  steady- 
state  solution  {ip,C)  was  obtained  by  DNS  and  the  time  i  at  which  the  damping  rate  of  residuals 
converged  to  within  a  predefined  tolerance  of  relative  discrepancy  10“®  between  successive  values 
of  a  and  the  results  for  ip  were  calculated.  The  value  of  a  as  well  as  the  relative  discrepancy 
Alp  =  \{ip{i)  —  ip)/ip\  between  the  estimated  and  the  exact  steady-states  is  also  shown;  the  level 
at  which  the  eigenmode  being  damped  is  present  in  the  transient  solution  at  time  t  may  be 
inferred  from  Aip. 

The  most  significant  result  of  this  table  is  the  ratio  iji.  The  case  Re  =  100  is  typical  of  one 
in  which  the  least-stable  eigenmode  determines  the  transient  behaviour  of  the  DNS  throughout 
most  of  the  time-integration  process.  With  the  results  for  a  converging  quite  quickly,  the  desired 
converged  steady-state  may  be  obtained  at  a  time  between  a  quarter  at  the  coarsest  and  a  fifth 
at  the  finest  resolution  of  the  time  required  by  the  time-marching  algorithm  for  the  residuals 
to  be  eliminated.  The  result  for  a  is  only  marginally  affected  by  resolution  and  time-step;  the 
precise  times  at  which  a  converges  are  affected  by  a  small  amount  when  reflning  the  grid,  with 
the  finest  resolution  results  converging  earlier.  In  all  cases  use  of  the  algorithm  of  §  6  (b )  v 
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results  in  substantial  savings  compared  with  the  otherwise  necessary  computing  effort.  The 
spatial  distribution  of  the  difference  Aip,  obtained  using  48  collocation  points  to  discretise  each 
spatial  direction,  is  shown  in  Fig.  36;  aside  from  the  level  of  Aip  it  is  interesting  to  notice  that 
the  discrepancy  between  the  two  solutions  attains  its  maximum  values  in  the  centre  of  the  cavity 
and  neither  the  singularity  of  the  boundary  conditions  nor  the  corner  vortices  are  manifested 
in  this  quantity.  The  same  qualitative  behaviour  was  shown  by  all  distributions  of  Aip  at  lower 
resolutions.  An  estimate  of  the  converged  solution  ip  obtained  by  application  of  (6.15)  at  f  =  15 
may  be  found  in  Fig.  37,  drawn  as  contours  at  the  levels  presented  by  Ghia  et  al.  (1982).  No 
cosmetic  post-processing  of  the  results  has  been  applied,  with  values  presented  at  the  collocation 
points  used.  As  it  is  to  be  expected  by  the  results  of  Table  8  the  agreement  between  ip  and  the 
result  of  Ghia  et  al.  (1982)  is  remarkable. 

At  Re  =  1000,  a  converges  at  approximately  the  same  fraction  of  total  integration  time  as  in 
the  Re  =  100  results.  However,  compared  with  the  Re  =  100  case  where  the  discrepancy  between 
estimated  and  converged  steady-states  is  three  to  four  orders  of  magnitude  smaller  compared 
with  that  shown  by  Aip,  here  only  one  order  of  magnitude  difference  between  the  maxima  of  Aip 
and  Alp  is  shown.  Though  small,  the  discrepancy  between  ip  and  ip  is  much  larger  than  roundoff 
level,  implying  that  elimination  of  the  least  stable  eigenmode  from  the  time-dependent  signal  for 
Ip  at  Re  =  1000  does  not  suffice  to  deliver  the  desired  ip.  Another  observation  that  may  be  made 
by  comparing  the  results  of  the  Re  =  100  and  Re  =  1000  cases  is  that  at  approximately  the  same 
value  of  iji  «  0.23,  Aip  is  higher  by  about  an  order  of  magnitude  at  Re  =  1000  compared  with 
that  at  Re  =  100.  In  searching  for  an  explanation  of  this  behaviour,  three  factors  may  be  recalled. 
First,  convergence  of  a  between  successive  values  is  a  necessary  but  not  sufficient  condition  for 
the  algorithm  of  §  6  (b )  v  to  deliver  accurate  results;  the  converged  a  should  be  compared 
with  the  corresponding  damping  rates  in  the  least-stable  part  of  the  eigenspectrum  of  the 
converged  steady-state  Q^.  Second,  as  the  Reynolds  number  increases  the  damping  rates  of  all 
global  eigenmodes  decrease,  suggesting  that  increasingly  longer  integration  times  are  necessary 
in  the  case  of  a  higher  Reynolds  number  in  order  for  the  residuals  to  subside  to  the  same  level  as 
in  a  lower  Reynolds  number  case.  Third,  the  separation  of  the  eigenvalues  in  the  global  spectrum 
plays  a  significant  role  in  attracting  the  transient  solution.  A  distinction  must  be  made  between 
the  early  and  the  late  stages  of  the  transient  behaviour  of  the  DNS  solution.  In  the  latter  it  is 
the  least-damped  eigenmode  which  must  eventually  be  damped  in  order  for  a  steady-state  to 
be  obtained.  During  the  early  stages  of  the  simulation,  on  the  other  hand,  an  arbitrary  initial 
condition  may  need  a  large  number  of  damped  global  eigenmodes  in  order  to  be  reconstructed.  It 
is,  therefore,  conceivable  that  at  the  early  stages  of  the  simulation  a  number  of  eigenmodes  other 
than  the  least-damped  one  are  present  in  the  transient  solution.  However,  as  time  progresses, 
increasingly  more  of  these  additional  eigenmodes  subside  on  account  of  their  large  damping 
rates,  to  the  effect  that  only  the  least-damped  mode  remains  to  determine  the  behaviour  of  the 
residual.  In  other  words,  as  time  progresses,  equation  (6.17)  reduces  to  (6.11)  and  the  theory  of 
§  6  (b )  V  focussing  on  a  single  damped  eigenmode  is  applicable. 

This  conjecture  may  easily  be  put  to  test  by  simply  permitting  the  time  integration  in  the 
second  DNS  to  proceed  beyond  t  while  monitoring  on  the  one  hand  a  against  Hi  and  on  the 
other  hand  Aip  in  the  process;  the  results  may  be  found  in  Table  9.  At  both  the  lower  and 
the  higher  Reynolds  number  further  integration  of  the  equations  of  motion  in  time  results  in 
all  but  the  least-stable  eigenmode  being  eliminated  from  the  signal,  as  clearly  demonstrated  by 
the  progressive  agreement  between  the  damping  rate  of  residuals  a  and  the  damping  rate  Hi  of 
the  least  stable  (A  =  0)— global  flow  eigenmode.  Consistent  with  this  result  is  the  increasingly 
improved  accuracy  by  which  the  algorithm  of  §  6  (b )  v  returns  the  estimate  of  the  converged 
steady  state,  as  shown  by  the  minimum  and  maximum  values  of  Aip  also  cited.  Interestingly,  ip 
may  be  recovered  at  the  same  low  level  of  discrepancy  in  the  two  Reynolds  number  cases,  f.e. 
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0(10“®)  at  Re  =  100,  i/i  =  .25  and  Re  =  1000,  i/i  =  .35,  although  the  agreement  of  a  with  Q, 
in  the  Re  =  100  is  about  an  order  of  magnitude  better  than  that  in  the  Re  =  1000  case. 

(v)  Three-dimensionality  as  a  consequence  oe  amplieied  (A  /  0)  two-dimensional 

LINEAR  EIGENMODES 

Finally,  we  turn  our  attention  to  the  differences  between  two-  and  three-dimensional  sim¬ 
ulations  on  account  of  growing  global  linear  instability  modes.  While  the  physics  behind  the 
instability  mechanisms  is  universal,  the  lid-driven  cavity  flow  example  serves  again  as  a  demon¬ 
strator,  with  the  differences  between  two-  and  three-dimensional  numerical  simulation  results  in 
this  flow  being  well  established  (Ku  et  al.  1987).  Here  we  call  upon  the  global  linear  instability 
theory  to  discuss  their  straightforward  explanation. 

It  is  possible  that  while  the  (A  =  0)— eigenmodes  at  a  certain  Reynolds  number  are  damped 
there  exist  unstable  A  /  0  global  flow  eigenmodes.  Indeed,  Ding  and  Kawahara  (1998)  have 
shown  that  at  Re  =  950  the  flow  is  unstable  to  modes  having  A  G  [A;,  A/j]  with  A;  =  27r/L/i  «  6.6 
and  A/i  =  27r/Lj  «  8.3,  while  the  domain  of  unstable  wavenumbers  systematically  broadens  in 
both  directions  on  the  A— axis  as  the  Reynolds  number  increases.  There  exist  two  possibilities  of 
introduction  of  three-dimensionality  by  means  of  DNS,  either  by  considering  spanwise  periodicity 
(pDNS)  or  by  taking  an  aperiodic  spanwise  domain  bounded  by  solid  walls  (aDNS).  In  the  case 
of  pDNS  the  integration  domain  in  the  spanwise  direction  is  defined  through  discrete  integer 
multiples  of  a  fundamental  wavenumber  Aq  such  that  ^2  =  271/ (nAo),  n  =  1, 2,  •  •  •,  while  in  aDNS 
Lz  is  a  continuous  free  parameter.  We  discuss  the  two  possibilities  separately;  in  both  cases  we 
restrict  the  discussion  to  simulations  performed  under  initial  and  boundary  conditions  such  that 
linear  instability  mechanisms  alone  can  drive  nonlinearity. 

If  a  three-dimensional  pDNS  is  performed  at  Re  =  950  and  a  spanwise  length  of  the  integration 
domain  is  chosen  such  that  Aq  >  8.3,  that  is  <  0.76,  neither  Aq  nor  any  of  the  harmonics 
of  this  global  linear  eigenmode  can  be  amplified.  As  a  consequence  one  may  predict,  without 
performing  the  three-dimensional  simulation,  that  the  latter  will  converge  in  time  to  the  same 
steady-state  solution  to  which  a  two-dimensional  {d/dz  =  0)  simulation  converges.  At  the  same 
Reynolds  number  value,  a  choice  of  spanwise  wavelength  G  [Li,  L/j]  =  [0.76, 0.95]  will  result  in 
exponential  amplification  and,  eventually,  turbulent  flow  on  account  of  the  unstable  fundamental 
wavenumber  which  is  implicitly  defined  by  a  spanwise  wavelength  within  this  range.  Finally,  if 
Lz  >  0.95  two  distinct  situations  may  be  obtained;  with  the  fundamental  wavenumber  being 
stable  (Aq  <  6.6),  may  be  taken  such  that  none  of  its  harmonics  fit  within  the  domain 
of  unstable  wavenumbers  at  this  Reynolds  number,  or  an  may  be  chosen  such  that  some 
harmonic  may  be  amplified.  While  in  the  first  case  the  two-dimensional  steady-state  solution 
will  be  obtained,  the  result  of  a  three-dimensional  simulation  in  the  second  case  will  be  transition 
to  a  turbulent  flow  state.  The  case  of  an  aDNS  may  be  perceived  as  a  special  case  of  a  pDNS, 
since  the  homogeneous  Dirichlet  conditions  imposed  on  the  disturbance  quantities  are  a  subset 
of  those  admissible  in  a  periodic  simulation.  Here  there  exist  two  possibilities,  depending  on 
whether  is  smaller  or  larger  than  Lj.  In  the  first  case  two-  and  three-dimensional  simulations 
will  deliver  identical  converged  steady-state  solutions  while  in  the  second,  which  includes  the 
well-studied  case  of  a  cubic  cavity,  transition  to  turbulence  should  be  expected  on  account  of  at 
least  one  three-dimensional  (A  /  0)  eigenmode  having  a  wavenumber  which  fits  into  A  G  [6.6, 8.3] 
at  this  Reynolds  number  value.  At  higher  Reynolds  number  values  the  situation  is  qualitatively 
analogous  for  aDNS,  with  the  dichotomy  in  wavenumbers  being  determined  by  the  highest 
neutrally  stable  wavenumber  value.  For  pDNS,  on  the  other  hand,  the  analogous  discussion  to 
that  at  Re  =  950  applies  at  Lz  <  Li  and  G  [Li,Lh\.  There  exists  a  Reynolds  number  value, 
though,  at  which  A/j  >  2Aj;  in  such  a  situation,  if  Lz  >  Lh  there  will  always  be  some  harmonic 
of  Aq  which  will  correspond  to  an  unstable  mode  having  A  =  nAo  G  [^D^/i]  which  will  be  liable 
to  linear  amplification  in  the  three-dimensional  simulation  and  eventual  departure  of  the  three- 
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from  the  two-dimensional  numerical  simulation  results.  The  lid-driven  cavity  with  its  large  body 
of  numerical  results  is  but  one  example  of  demonstration  of  this  behaviour;  we  are  currently 
applying  these  ideas  in  the  compressible  elliptic  cone  flowfield. 

{d)  Summary  and  Epilogue 

The  questions  raised  in  §  6  (a  )  may  be  answered  within  the  unifying  framework  of  global  linear 
instability  analysis  of  a  two-dimensional  steady  solution  of  the  equations  of  motion.  Aided  by  the 
results  of  a  numerically  well-studied  incompressible  flow  problem  we  were  able  to  attach  physical 
significance  to  the  transient  behaviour  of  two-dimensional  time-dependent  incompressible  direct 
numerical  simulation  results.  What  is  commonly  known  as  residual  in  the  simulation  is  either 
the  least  damped  two-dimensional  (A  =  0)— linear  eigenmode  of  the  converged  steady  state 
itself,  or  can  be  related  to  a  small  number  of  the  least  damped  modes  of  the  full  eigenvalue 
spectrum.  As  the  Reynolds  number  increases,  all  global  two-dimensional  eigenmodes  become 
increasingly  less  damped,  until  a  parameter  value  is  reached  beyond  which  no  steady-state 
solution  exists.  The  physical  information  which  is  suppressed  in  two-dimensional  simulations 
based  on  the  steady  formulation  of  the  equations  of  motion  concerns  the  dynamical  behaviour 
of  these  two-dimensional  linear  eigenmodes.  While  unsteadiness  should  not  be  interpreted  as 
amplification  of  the  global  linear  (A  =  0)— eigenmodes,  on  the  simple  grounds  of  the  absence  of  a 
converged  steady-state  upon  which  the  latter  would  develop,  the  process  leading  to  unsteadiness 
is  directly  linked  with  the  diminishing  magnitude  of  damping  rates  of  the  global  linear  modes 
as  the  flow  Reynolds  number  increases,  and  the  associated  prevalence  of  nonlinearity. 

When  a  steady-state  solution  exists,  the  insight  gained  from  the  association  of  the  transient 
behaviour  in  two-dimensional  DNS  with  the  results  of  the  nonparallel  linear  instability  anal¬ 
ysis  of  the  converged  steady-state  may  be  utilised  in  a  threefold  manner.  First,  an  algorithm 
may  be  constructed,  to  recover  the  steady-state  solution  from  transient  data  taken  well  before 
convergence,  thus  making  further  time- integration  of  the  equations  of  motion  redundant.  The 
algorithm,  whose  building  elements  were  presented  in  §  6  (b )  v,  is  based  on  identification  of 
the  parameters  pertaining  to  the  linear  eigenmodes  which  determine  the  transient  behaviour  of 
the  solution,  namely  the  damping  rate  a  and  the  frequency  Dr  of  the  least  stable  eigenmodes. 
Results  shown  in  §  6  (c )  iv  on  the  example  problem  studied  have  demonstrated  that  up  to 
three-quarters  of  the  otherwise  necessary  computing  effort  may  be  saved  by  application  of  the 
theory  of  §  6  (b )  v.  Second,  the  results  of  a  nonparallel  linear  instability  analysis  of  the  con¬ 
verged  steady-state  can  be  used  as  a  quality  test  of  the  obtained  solution,  if  the  latter  has  been 
obtained  using  a  time-accurate  solution  approach.  The  rate  of  decay  of  the  residual  which  ulti¬ 
mately  has  to  be  damped  in  order  for  a  converged  steady-state  to  be  obtained  should  equal  the 
damping  rate  of  the  least-stable  eigenmode,  if  both  numbers  are  substantially  larger  than  zero 
in  magnitude.  Disagreement  of  these  two  quantities  indicates  that  the  obtained  steady-state  still 
contains  an  unsteady  component  which  must  be  eliminated  by  further  time-integration,  or  by 
application  of  the  ideas  of  §  6  (b )  v.  Third,  the  time  necessary  for  the  reduction  of  residuals  to 
machine-roundoff  level  may  also  be  estimated  using  nonparallel  linear  instability  theory  and  is 
inversely  proportional  to  the  damping  rate  of  the  least  damped  linear  eigenmode.  Using  the  value 
of  the  damping  rate  obtained  by  extrapolation  of  data  at  lower  Reynolds  numbers  one  predicts 
that  in  the  square  lid-driven  cavity  at  Re  =  10‘*  a  steady-state  solution,  if  one  exists,  may  be 
obtained  after  integrating  the  unsteady  equations  of  motion  for  time  in  excess  of  t  =  4000  as 
calculated  from  (6.10)  and  non-dimensionalised  with  the  lid- velocity  and  the  cavity  length. 

Well  before  the  flow  tends  to  lose  its  stability  to  two-dimensional  linear  eigenmodes,  three- 
dimensional  (A  /  0)— disturbances  may  be  amplified.  Depending  on  the  size  of  the  observa¬ 
tion  window  in  the  third  spatial  dimension,  this  amplification  of  three-dimensional  global  dis¬ 
turbances  can  explain  the  differences  between  the  results  of  the  two-  and  three-dimensional 
DNS.  Again,  caution  is  warranted  at  this  point  not  to  confuse  amplification  of  the  global,  two- 
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dimensional  instabilities  discussed  here  with  solutions  of  the  classic  ordinary-differential-equation 
based  eigenvalue  problem,  which  are  incorporated  in  those  of  (6. 6-6. 9);  both  mechanisms  may 
provide  amplification,  as  the  laminar  separation  bubble  flow  example  has  clearly  demonstrated 
(Theofilis  et  al.  2000).  Conversely,  nonparallel  linear  instability  theory  provides  a  handle  to 
probe  into  the  physics  of  the  flow  in  (three-dimensional)  physical  space  using  two-dimensional 
DNS  results,  before  resorting  to  computationally  intensive  three-dimensional  spatial  DNS,  at 
least  as  far  as  the  response  of  the  flow  to  small-amplitude  excitations  is  concerned.  Solution  of 
the  partial-derivative  eigenvalue  problem  not  only  answers  the  question  whether  new  physics  is  to 
be  learnt  by  performing  the  three-dimensional  DNS  at  a  given  set  of  parameters  but  also  provides 
information  on  the  physical  mechanism  which  leads  flow  to  deviate  from  two-dimensionality. 

Based  on  the  findings  presented  we  may  extend  the  discussion,  in  the  form  of  proposed  future 
work,  to  both  one  and  three  nonperiodic  spatial  directions.  Both  an  one-dimensional  and  a 
three-dimensional  steady-state  solution  Qj  may  be  recovered  by  application  of  the  ideas  discussed 
herein  for  the  case  of  two  nonperiodic  spatial  directions.  In  the  case  of  an  one-dimensional  profile 
Qb  being  sought  by  time-marching  the  equations  of  motion,  taking  two  spatial  directions  as 
periodic  and  resolving  the  third,  the  associated  linear  instability  problem  to  be  solved  is  based  on 
the  classic  system  of  the  one-dimensional  Orr-Sommerfeld  and  Squire  linear  instability  equations 
to  which  (6. 6-6. 9)  reduce  if  the  dependence  of  the  basic  flow  on  one  of  the  two  resolved  spatial 
directions,  say  x,  is  neglected  such  that  this  spatial  direction  may  be  taken  as  homogeneous  as 
far  as  the  disturbance  field  is  concerned.  The  linear  mode  associated  with  the  residuals  is  the 
least  stable  member  of  the  spectrum  obtained  at  k  =  A  =  0,  k  and  A  being  the  wavenumbers 
along  the  periodic  spatial  directions,  x  and  z.  It  is  well  appreciated  in  this  case  that  agreement  of 
the  time-accurate  simulation  results  and  those  of  the  one-dimensional  linear  instability  problem 
is  a  minimum  simulation  quality  criterion  (Kleiser  and  Schumann  1980;  Canuto  et  al.  1987). 
However,  given  current  hardware  capabilities,  it  is  likely  that  an  one-dimensional  Qj  will  be 
sought  by  a  direct  algorithm,  rather  than  by  time-marching  the  unsteady  equations  of  motion. 

An  extension  of  the  algorithm  presented  for  the  recovery  of  a  two-dimensional  is  also  possi¬ 
ble  in  the  case  of  flow  developing  in  three  nonperiodic  spatial  directions.  In  this  case  the  existence 
of  a  steady-state  Qj  is  synonymous  with  stability  of  all  eigenmodes  of  the  flow  but  current  hard¬ 
ware  technology  makes  the  solution  of  the  corresponding  three-dimensional  partial  derivative 
eigenvalue  problem  impractical.  On  the  other  hand  the  ideas  presented  in  section  6  (b )  v  may 
be  used  in  order  to  recover  a  three-dimensional  steady  state  once  a  regime  of  linear  damping 
of  residuals  has  been  identified.  This  is  of  immediate  interest  in  the  case  of  the  elliptic  cone 
geometry,  where  recovery  of  the  converged  steady  laminar  three-dimensional  steady  state  by 
application  of  the  ideas  presented  in  this  section  is  currently  being  actively  pursued. 
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7.  Discussion 

This  report  commences  presentation  of  our  efforts  towards  global  linear  instability  analysis  of 
high-speed  flow  around  elliptic  cones.  We  have  presented  the  framework  within  which  the  analysis 
is  envisaged  to  be  performed  and  it  has  been  argued  that  a  global  linear  analysis  in  which  all 
three  spatial  directions  are  resolved  is  most  efficiently  performed  by  DNS.  Our  ultimate  concern 
is  with  a  global  linear  analysis  in  which  two  spatial  directions  are  resolved  while  the  third  is 
taken  to  be  homogeneous.  Strictly,  no  such  homogeneous  spatial  direction  exists  in  the  flow  over 
an  elliptic  cone;  however  the  dependence  of  the  basic  flow  quantities  on  the  directions  normal  to 
the  cone  surface  dominates  over  that  along  the  cone  generators  and  the  global  analysis  may  be 
performed  by  resolving  the  former  two  and  neglecting  the  third  spatial  direction  when  addressing 
the  eigenvalue  problem.  Instead  of  making  this  simplifying  assumption  when  computing  the  basic 
flow  also,  we  have  expanded  the  scope  of  our  original  intention  to  perform  a  two-dimensional 
solution  for  the  basic  flow  and  have  opted  for  a  numerical  solution  of  the  full  three-dimensional 
problem.  We  chose  the  AUSMDV  shock-capturing  scheme  for  the  discretisation  of  the  convective 
fluxes  in  the  equations  of  motion.  The  scheme  was  validated  and  its  satisfactory  performance  in 
resolving  flow  discontinuities  was  established  before  it  was  to  the  elliptic  cone  geometry. 

Solutions  of  the  steady  three-dimensional  equations  of  motion  were  obtained  next.  A  half- 
cone  model  was  considered  and  symmetry  was  imposed  on  the  solution,  implying  that  all  the 
results  obtained  herein  pertain  to  angle  of  bank  /3  =  0.  The  elliptic  cone  was  terminated  by  a 
planar  surface  and  the  object  was  embedded  into  three-dimensional  space  which  extended  well 
away  from  the  surface  of  the  elliptic  cone  object.  Terminating  the  cone  with  a  planar  surface 
resulted,  as  expected,  in  large  flow  separation  at  the  base  both  in  subsonic  and  supersonic  flow. 
In  forthcoming  basic  flow  calculations  we  intend  to  investigate  a  compound  object  composed  of 
the  same  elliptic  cone  and  a  half  prolate  spheroid/ellipsoid  joined  together  so  as  to  minimise  the 
curvature  jump  at  base  and  thus  reduce  the  intensity  of  separation!  behind  the  cone. 

Instead  of  devoting  our  efforts  entirely  to  the  refinement  of  the  solution  at  a  single  set  of 
parameters,  the  approach  taken  in  generating  solutions  of  the  equations  of  motion  was  the 
creation  of  a  database  of  initial  conditions  at  different  parameters.  This  approach  is  expected  to 
assist  subsequent  global  linear  instability  studies  by  providing  appropriate  attractors  for  fine- 
resolution  basic  flow  calculations,  as  required.  In  order  to  ensure  laminar  flow  we  have  kept  a 
constant  low  Reynolds  number  value  Re  =  10^  and  permitted  variation  of  the  angle  of  attack 
and  the  flow  Mach  number.  Subsonic  flow  solutions  were  obtained  at  M  =  0.5,  a  =  10°  and 
M  =  0.5,  a  =  20°  while  the  same  hybrid  grid  permitted  recovery  of  supersonic  solutions  at 
M  =  2,  a  =  20°.  At  M  =  4,  a  =  20°  the  sharpness  of  the  gradients  developing  called  for  an 
adaptation  of  the  calculation  grid,  firstly  using  the  solution  gradients  to  redistribute  the  available 
gridpoints  in  three-dimensional  space  and  subsequently  increasing  the  number  of  points  by  some 
30%  compared  with  those  at  the  lower  Mach  number  values.  The  same  trend  has  been  observed 
at  preliminary  studies  at  M  >  4,  where  even  higher  resolutions  to  those  used  herein  are  expected 
to  be  necessary. 

We  then  turned  our  attention  to  the  instability  analysis  of  the  recovered  flowfields.  It  has  been 
argued  that  a  global  three-dimensional  linear  instability  analysis,  in  which  a  converged  steady 
three-dimensional  basic  flow  forms  the  variable  coefficients  of  the  partial-differential-equation 
three-dimensional  eigenvalue  problem,  is  uninteresting  from  a  physical  point  of  view  in  its  own 
right.  Instead,  we  have  focussed  our  attention  on  the  two-dimensional  eigenvalue  problem  and 
presented  the  equations  governing  alternative  forms  of  inviscid  global  two-dimensional  linear 
analysis.  In  the  limit  of  a  =  0  we  have  derived  the  extension  of  the  Rayleigh  equation  governing 
global  inviscid  instability  of  compressible  flows.  In  case  the  analysis  of  a  flowfield  in  which  a  /  0 


t  itself  a  source  of  global  instability  (Theofilis  et  al.  2000) 
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is  of  interest,  it  has  been  argued  that  no  numerical  advantage  exists  in  considering  an  inviscid 
analysis;  the  global  instability  studies  is  then  best  performed  in  a  viscous  framework. 

Finally,  in  view  of  the  intimate  connection  between  global  instability  analysis  and  steady- 
state  solutions  of  the  equations  of  motion,  we  have  devoted  substantial  efforts  in  identifying 
residuals  in  time-accurate  simulations  as  the  least  stable  global  eigenmodes  of  the  converged 
in  time  steady-state  solution.  This  knowledge,  independently  verified  on  the  classic  lid-driven 
cavity,  permits  utilising  global  linear  theory  to  construct  a  theoretically-founded  convergence 
acceleration  technique  towards  a  steady-state  solution  of  the  equations  of  fluid  flow  motion.  An 
algorithm  has  been  presented  using  which  the  desired  steady  state  may  be  recovered  by  simple 
algebraic  operations  on  transient  data  taken  well  before  the  time-integration  procedure  has 
converged.  In  the  worked  example  substantial  savings  in  computing  effort  have  been  materialised 
compared  with  the  otherwise  necessary  time-integration  until  ’residuals’  are  reduced  to  machine 
roundoff  level.  Extension  of  this  idea  on  a  three-dimensional  basic  flow  is  currently  underway 
using  the  elliptic  cone  geometry.  Results  will  be  reported  in  due  course. 
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Variable 

Left  state  (subscript  l) 
0.0  <  :e  <  0.1 

Middle  state  (subscript  m) 
0.1  <  :e  <  0.9 

Right  state  (subscript  jj) 

0.9  <  :e  <  1.0 

P 

1 

1 

1 

u 

0 

0 

0 

P 

10® 

1 

o 

1—1 

10® 

Table  1.  Initial  conditions  for  the  blast-wave  problem 
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A.X  X  10'*  pup 


10 

0.1624 

X 

0.0116 

=  0 

73.9501 

0.0001 

5 

0.1580 

0.0057 

72.4378 

0.0000 

2 

0.1531 

0.0023 

70.3934 

0.0000 

1 

0.1507 

0.0012 

69.3482 

0.0000 

0.5 

0.1493 

0.0006 

68.8574 

0.0000 

10 

0.3201 

X 

-0.0003 

=  1 

19.2320 

0.0000 

5 

0.3170 

-0.0002 

19.2811 

0.0000 

2 

0.3140 

-0.0001 

19.2941 

0.0000 

1 

0.3124 

-0.0000 

19.2941 

0.0000 

0.5 

0.3113 

-0.0000 

19.2939 

0.0000 

Table  2.  Convergence  of  the  wall  values  of  the  primitive  variables  in  the  blast  wave  problem 
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Re  =  400 

Ghia  et  al.  (1982) 

Schreiber  and  Keller  (1983) 

present  results 

Primary 

iP 

c 

{x,y) 

-0.1139 

2.29469 

(0.5547,  0.6055) 

-0.1140 

2.281 

(0.5571,0.6071) 

-0.1139 

2.29584 

(0.5535,0.6054) 

LL 

ip 

C 

{x,y) 

1.42  X  10“® 
-0.0570 
(0.0508,  0.0469) 

1.45  X  10“® 

-0.0471 

(0.0500,0.0429) 

1.40  X  10“® 
-0.05685 
(0.0510,0.0466) 

LR 

ip 

C 

{x,y) 

6.42  X  10“^ 
-0.4335 
(0.8906,0.1250) 

6.44  X  10“^ 

-0.394 

(0.8857,0.1143) 

6.41  X  10“^ 
-0.44802 
(0.8852,0.1217) 

Re  =  1000 

Ghia  et  al.  (1982) 

Schreiber  and  Keller  (1983) 

present  results 

Primary 

ip 

C 

(®,2/) 

-0.117929 

2.04968 

(0.5313,0.5625) 

-0.11603 

2.02600 

(0.52857,0.56429) 

-0.118902 

2.068251 

(0.529654,0.565018) 

LL 

iP 

c 

{x,y) 

2.31  X  10“^ 
-0.36175 
(0.0859,  0.0781) 

2.17  X  10“® 

-0.302 

(0.08571,0.07143) 

2.354097  X  10“^ 
-0.337187 
(0.081549,0.077839) 

LR 

ip 

C 

{x,y) 

1.75  X  10“® 
-1.15465 
(0.8594,  0.1094) 

1.70  X  10“® 

-0.999 

(0.86429,0.10714) 

1.744028  X  10“® 
-1.097921 
(0.867381,0.114469) 

Re  =  3200 

Primary 

UL 

LL 

LR 

GGS 

ip 

C 

{x,y) 

-0.12038 

1.98860 

(0.5165,  0.5469) 

7.27682  X  10“^ 
-1.71161 
(0.0547,  0.8984) 

9.7823  X  10“^ 
-1.06301 
(0.0859,0.1094) 

3.14  X  10“® 
-2.27365 
(0.81255,0.0859) 

present 

ip 

C 

{x,y) 

-0.12181 

1.961154 

(0.51722,0.54089) 

7.11201  X  10“^ 
-1.65335 
(0.0524,0.8981) 

1.12331  X  10“® 
-1.16397 
(0.08106,0.12052) 

2.82648  X  10“® 
-2.24381 

(0.82281,0.084648) 

Re  =  4000 

Primary 

UL 

LL 

LR 

SK 

ip 

C 

{x,y) 

-0.11237 

1.805 

(0.51875,  0.53750) 

1.12  X  10“® 
-1.067 

(0.08125,0.11875) 

2.80  X  10“® 
-2.145 

(0.81875,0.07500) 

present 

ip 

C 

{x,y) 

-0.12203 

1.94949 

(0.51597,  0.53846) 

1.073  X  10“® 

-1.91234 

(0.06098,0.90387) 

1.24736  X  10“* 
-1.27899 
(0.08055,0.12482) 

2.95426  X  10“* 
-2.42032 
(0.81640,0.07983) 
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Table  3.  Location  and  value  of  the  maxima  of  the  primary  and  the  lower-left  (LLJ,  lower-right  (LR)  and  upper-left 
(UL)  secondary  vortices  in  the  steady  state  solution  for  ^  and  ^  at  Re  =  400, 1000,  3200  and  4000;  comparisons 
with  Ghia  et  al.  (1982)  (GGS)  and  Schreiber  and  Keller  (1983)  (SK). 
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Primary 

i’ 

Re  = 

UL 

c 

:  400 

LL 

loV  C 

LR 

loV  C 

GGS 

-0.113909 

2.29469 

1.41951 

-0.05697 

6.42352 

-0.4335 

present 

-0.113989 

2.29463 

1.47210 

-0.05711 

6.42406 

-0.4329 

SK 

-0.11399* 

2.2898* 

1.45 

-0.04710 

6.440 

-0.3940 

present 

-0.113982 

2.29184 

1.40 

-0.04766 

6.373 

-0.4030 

Re  = 

1000 

Primary 

UL 

LL 

LR 

C 

i’ 

c 

loV 

c 

lOV 

c 

GGS 

-0.117929 

2.04968 

2.31129 

-0.36175 

1.75102 

-1.1547 

present 

-0.118902 

2.06839 

2.37806 

-0.36575 

1.77911 

-1.1486 

SK 

-0.11894* 

2.0677* 

2.1700 

-0.302000 

1.700 

-0.9990 

present 

-0.118905 

2.068234 

2.3151 

-0.312162 

1.763 

-1.0481 

Re  = 

3200 

Primary 

UL 

LL 

LR 

c 

loV 

c 

lO^ip 

c 

lO^ip 

C 

GGS 

-0.120377 

1.9886 

7.27 

-1.71161 

0.98 

-1.06301 

-2.27365 

present 

-0.121777 

1.9612 

7.08 

-1.73137 

1.09 

-1.00607 

-2.25511 

Re  = 

4000 

Primary 

UL 

LL 

LR 

c 

i’ 

c 

lOV 

c 

lOV 

c 

SK 

-0.12202* 

1.9498* 

1.1200 

-1.0670 

2.8000 

-2.14500 

present 

-0.122026 

1.94960 

1.2411 

-1.1427 

2.9228 

-2.31944 

Table  4.  Comparison  of  the  interpolated  values  of  our  solutions  on  the  maxima  presented  by  Ghia  et 
and  Schreiher  and  Keller  (1983).  An  asterisk  denotes  Richards  on- extrapolated  data  in  the  latter 


al.  (1982) 
work. 
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RH  DK  present  results 


A 

Qr 

Qi 

Qi 

fir 

fli 

1 

±0.00 

-0.34 

-0.3183 

±0.0000 

-0.3297 

2 

±0.00 

-0.23 

-0.2248 

±0.0000 

-0.2267 

3 

±0.11 

-0.29 

-0.2924 

±0.1073 

-0.2954 

4 

±0.28 

-0.30 

-0.2969 

±0.2810 

-0.2956 

5 

±0.43 

-0.34 

-0.3431 

±0.4260 

-0.3404 

6 

±0.58 

-0.39 

-0.3893 

±0.5821 

-0.3844 

7 

±0.67 

-0.41 

-0.4073 

±0.6733 

-0.4013 

8 

±0.72 

-0.45 

-0.4637 

±0.7232 

-0.4587 

9 

±0.76 

-0.54 

-0.5504 

±0.7622 

-0.5473 

Table  5.  Comparison  of  the  least  stable  eigenmode  at  Re  =  200  against  the  results  of  Ramanan  and  Homsy 
(1994)  (RH)  and  the  graphically  (digitally)  reproduced  growth  rate  result  of  Ding  and  Kawahara  (1998)  (DK). 
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Re 

100 

200 

300 

Resolution 

a 

<7 

a 

16  X  16 

-0.5404 

-0.3248 

-0.2865 

24  X  24 

-0.5407 

-0.3319 

-0.2843 

32  X  32 

-0.5409 

-0.3318 

-0.2842 

40  X  40 

-0.5409 

-0.3318 

-0.2842 

Qi 

-0.5410 

-0.3319 

-0.2845 

Table  6.  Numerical  results  for  the  rate  of  decay  of  residuals  a  as  a  function  of  resolution  at  different  low  Reynolds 
numbers.  Also  shown  the  result  of  numerical  solution  of  (6. 6-6. 9)  for  the  imaginary  part  of  the  eigenvalue  Oi, 
using  the  respective  converged  steady-state  as  basic  flow. 
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Re 

Qi 

X  100 

<T 

2500 

-0.0253 

1.2 

5000 

-0.0112 

8.9 

7500 

-0.0093 

17.8 

Table  7.  Numerical  results  for  the  damping  rate  Oi  of  the  least  stable  (A  =  0)  — eigenmode  at  Re  =  2500,  5000  and 
7500  and  its  discrepancy  in  percentage  terms  from  the  rate  of  decay  of  residuals  a. 
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Re  =  100 

Resolution 

16  X  16 

24  X  24 

32  X  32 

At 

0.01 

0.01 

0.005 

t 

50.43 

50.42 

49.005 

i 

12.71 

12.79 

11.07 

—a 

0.540246 

0.540214 

0.540876 

max(Ai^) 

5.3(-8) 

8.8(-7) 

4.6(-6) 

min(Ai^) 

3.6(-9) 

7.9(-8) 

5.1(-7) 

max(A^) 

3.4(-4) 

3.5(-4) 

1.0(-3) 

min(A^) 

1.5(-5) 

1.2(-3) 

1.6(-2) 

Re  =  1000 

Resolution 

24  X  24 

32  X  32 

40  X  40 

At 

0.01 

0.01 

0.01 

t 

325.93 

323.61 

324.33 

i 

77.82 

77.53 

77.81 

— <7 

0.065808 

0.065657 

0.065336 

max(Ai^) 

2.9(-6) 

9.5(-5) 

3.1(-5) 

min(Ai^) 

4.9(-7) 

3.2(-6) 

3.3(-6) 

max(A^) 

3.7(-4) 

3.6(-4) 

2.5(-4) 

min(A^) 

3.7(-4) 

5.1(-4) 

3.3(-4) 

Table  8.  Recovery  of  ip  from  transient  data  at  Re  =  400  and  1000  as  function  of  resolution  and  time-step  used  in 

the  DNS.  x{y)  =  x  x  10®. 
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1  X  100 

max(A^) 

Re 

min(A^) 

=  100 

s/  1  fin 

22.59 

4.6(-6) 

5.1(-7) 

0.023 

25.02 

4.0(-8) 

1.9(-8) 

0.018 

30.90 

2.5(-9) 

1.2(-9) 

0.011 

35.12 

1.8(-10) 

I.O(-IO) 

0.009 

Re 

=  1000 

23.99 

3.1(-5) 

3.3(-6) 

3.97 

25.02 

1.3(-5) 

1.5(-6) 

3.35 

30.23 

2.1(-6) 

1.2(-7) 

0.65 

35.02 

1.5(-8) 

5.8(-9) 

0.24 

Table  9.  Recovery  o/Qj  at  several  Reynolds  numbers  from  transient  data  at  times  beyond  that  at  which  to  converges. 
The  discrepancy  between  a  and  Oi  of  the  least  stable  eigenmode  is  also  presented.  Re  =  100  run  on  a  32^  grid; 
Re  =  1000  run  on  a  40^  grid.  x{y)  =  x  x  10®. 
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Figure  2.  Global  view  of  the  grid  utilised.  Only  the  surface  and  farfield  discretisation  are  shown. 
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Figure  3.  Detail  of  the  near-field  discretisation  of  the  elliptic  cone  surface  and  the  symmetry  plane,  viewing 
downstream  the  cone.  The  hybrid  grid  can  be  seen  clearly  on  the  symmetry  plane. 
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Figure  4.  Detail  of  the  near-field  discretisation  of  the  elliptic  cone  surface  and  the  symmetry  plane,  viewing 
upstream  the  cone.  The  discretisation  of  the  cone  base  is  also  seen  in  this  plot. 
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AUSMDV  forRiemann’s  1st  Problem 


X 


AUSMDV  for  Riemann’s  2nd  Problem 


X 


Figure  5.  Two  shock-tube  problems  solved  by  an  AUSMDV  flux-vector  splitting.  The  primitive  variables  p,u,p 
and  T  are  presented  for  Sod’s  problem  in  the  upper  figure  and  for  Lax’s  problem  in  the  lower  figure  at  times  f  =  2 
and  1.3,  respectively.  In  both  cases  a  uniform  grid  with  Ax  =  1/100  has  been  used. 
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U 

175.659 

160.12 

144.58 

129.041 

113.501 

97.962 

82.4226 

66.8832 

51.3438 

35.8044 

20.265 

4.72564 

-10.8138 

-26.3532 


Figure  9.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  0.5,  a  =  10“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  a— wise  velocity  component  u  on  the  planes 
X  =  0.7  and  y  =  0- 
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FrameOOl  I  27  0ct2000  {  Grld:eRlpllc cone.grld, Pointdata:cone.pvaL10000 


V 

35.6133 

30.6584 

25.7034 

20.7484 

15.7934 

10.8385 

5.88348 

0.928506 

-4.02647 

-8.98144 


Figure  10.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  0.5,  a  =  10“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  2/— wise  velocity  component  v  on  the  plane 
X  =  0.7.  The  condition  n  =  0  is  satisfied  on  the  symmetry  plane  on  account  of  the  symmetry  imposed. 
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37.6873 

31.8522 

26.0171 

20.1819 

14.3468 

8.51163 

2.67649 

-3.15866 


Figure  11.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10^,  M  =  0.5,  a  =  10®;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  z— wise  velocity  component  w  on  the  planes 
X  =  0.7  and  y  =  0. 
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U 

170.289 

143.066 

115.843 

88.6194 

61.3963 

34.1731 

6.94999 

-20.2731 


U 


■ 


170.289 

143.066 

115.843 

88.6194 

61.3963 

34.1731 

6.94999 

-20.2731 


Figure  12.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  0.5,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  a— wise  velocity  component  u  on  the  planes 
X  =  0.7  and  y  =  0- 
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V 


61.0406 

48.9929 

36.9452 

24.8975 

12.8498 

0.802052 

-11.2456 

-23.2933 


Figure  13.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  0.5,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  2/— wise  velocity  component  v  on  the  plane 
X  =  0.7.  The  condition  n  =  0  is  satisfied  on  the  symmetry  plane  on  account  of  the  symmetry  imposed. 
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Figure  14.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  0.5,  a  =  20“ ;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  z— wise  velocity  component  w  on  the  planes 
X  =  0.7  and  y  =  0- 
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density 
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density 
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0.679217 


Figure  15.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  2,a  =  20°;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  p  on  the  planes  x  =  0.7  and  3/  =  0. 
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u 

■  614.739 

■  577.624 

■  540.509 

■  503.395 
=  466.28 
=  429.165 
—  392.051 
=  354.936 
H  317.821 
=  280.707 
=  243.592 
=  206.477 
=  169.362 
=  132.248 
H  95.1331 

■  58.0184 

■  20.9037 

U 

614.739 
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466.28 
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392.051 

354.936 
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132.248 

95.1331 

58.0184 

20.9037 


Figure  16.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  2,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  a— wise  velocity  component  u  on  the  planes 
X  =  0.7  and  y  =  0- 
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V 


158.098 

135.222 

112.345 

89.4687 

66.5923 

43.7158 

20.8393 

-2.03717 

-24.9137 

-47.7901 


Figure  17.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10*,M  =  2,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  j/— wise  velocity  component  v  on  the  plane 
X  =  0.7.  The  condition  n  =  0  is  satisfied  on  the  symmetry  plane  on  account  of  the  symmetry  imposed. 
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Figure  18.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  2,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  z— wise  velocity  component  w  on  the  planes 
X  =  0.7  and  y  =  0- 
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pressure 
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Figure  19.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  2,a  =  20“ ;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  p  on  the  planes  x  =  0.7  and  i/  =  0. 
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Figure  20.  A  qualitatively  self-similar  flowfield  pattern  of  the  z— wise  velocity  component  w  obtained  on  the 

elliptic  cone  at  M  =  2,a  =  20“. 
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2.92638 
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2.08439 
1 .6634 
1 .24241 
0.821419 
0.400427 


Figure  21.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10®,  M  =  A,a  =  20°;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  p  on  the  planes  x  =  0.7  and  y  =  0- 
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Figure  22.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10*,M  =  4,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  a— wise  velocity  component  u  on  the  planes 
X  =  0.7  and  y  =  0- 
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Frame  001  I  27  Oct  2000  I  Grld:eRlpllc cone.grld,  Pointdata:cone.pvaL40000 


V 

272.482 
216.581 
160.681 
^  104.781 
48.8809 
-7.0192 
-62.9193 
-118.819 


Figure  23.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10*,M  =  4,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  i/— wise  velocity  component  v  on  the  plane 
X  =  0.7.  The  condition  n  =  0  is  satisfied  on  the  symmetry  plane  on  account  of  the  symmetry  imposed. 
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247.458 

174.798 

102.137 

29.4767 

-43.1838 


Figure  24.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10*,M  =  4,  a  =  20“;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  the  z— wise  velocity  component  w  on  the  planes 
X  =  0.7  and  y  =  0. 
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Figure  25.  An  aspect-ratio  A  =  3  elliptic  cone  model  at  Re  =  10* ,  M  =  A,  a  =  20°;  shown  are  contour  lines 
equidistributed  between  the  minimum  and  maximum  values  of  p  on  the  planes  x  =  0.7  and  i/  =  0. 
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Figure  26.  Convergence  history  of  stream  function  ^(0.5,  0.5)  against  time  (left)  and  slope  of  this  curve  (right). 
Lower  to  upper  curves,  Re  =  100,  200  and  300,  respectively. 
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Figure  27.  Convergence  history  of  ^’(0.5,  0.5)  against  time  at  Re  =  500  (upper  left)  and  its  slope  (upper  right); 
lower  left  and  right,  respectively,  the  corresponding  results  at  Re  =  1000.  In  both  cases  superimposed  and  denoted 
by  symbols  are  the  eigenvalues  of  the  two  least  stable  stationary  modes. 
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Figure  28.  The  dependence  of  the  function  d(ln^*)/df  on  time  t,  showing  the  exponential  decay  of  a  single  travelling 
mode  (Or  «  0.97  ±  0.01)  superimposed  upon  the  least  damped  exponentially  decaying  stationary  disturbance  at 


Re  =  2500. 
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Figure  31.  The  correspondence  of  the  frequencies  of  the  damped  linear  (A  =  0)  two-dimensional  eigenmodes  of 
the  converged  steady-states  at  different  Reynolds  numbers  and  those  obtained  from  discrete  Fourier  transforms 
of  the  DNS  signals.  Re  =  2500. 
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Figure  35.  The  dependence  of  the  damping  rate  Oi  of  the  least  damped  two-dimensional  eigenmode  of  the  converged 
steady-state  at  a  Reynolds  number  on  Re  as  predicted  by  the  model  (6.21)  denoted  by  the  solid  line,  and  as 
calculated  by  numerical  solution  of  the  eigenvalue  problem  (6. 6-6. 9)  denoted  by  the  symbols. 


Contract  No.  P61775-99-WE049 


92 


V.  Theofilis 


0  0 


Figure  36.  The  spatial  distribution  of  the  difference  A-il>(x,y)  =  tj)  —  tj)  Re  =  100  using  Nx  =  Ny  =  48  Jacobi 

collocation  points. 
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Figure  37.  An  estimate  of  the  converged  solution  tj)  Re  =  100  obtained  by  evaluating  (6.15)  at  f  =  15  and 
using  Nx  =  Ny  =  48  Jacobi  collocation  points.  Iso-contours  are  drawn  at  the  levels  shown  by  Ghia  et  al.  (1982) 
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